ADRECMO

Method

This study is a monocentric, prospective, observational investigation conducted in Nancy Hospital between 10.10.2017 and 24.05.2022, enrolling 50 adult patients with refractory cardiogenic shock requiring veno-arterial extracorporeal membrane oxygenation (VA-ECMO) support (NCT03327493).

The study was carried out in accordance with the Declaration of Helsinki and approved by the local institutional ethics committee. Written informed consent was obtained from all patients or their legal representatives prior to inclusion.

Patient Population

Cardiogenic shock was defined as a systolic blood pressure (SBP) <90 mmHg or a mean arterial pressure (MAP) <65 mmHg despite adequate volume resuscitation, accompanied by signs of peripheral hypoperfusion (e.g., cold extremities, oliguria, altered mental status) and a cardiac index <2.2 L/min/m². Refractoriness was defined as a hyporesponsiveness to norepinephrine and/or the persistence of profound clinical signs of hypoperfusion despite optimal resuscitation efforts. Eligible patients were adults (≥18 years old) affiliated with a health insurance system. Cardiogenic shock etiologies included acute ischemic coronary artery disease, primary or ischemic-related dilated cardiomyopathy, myocarditis and stress-induced cardiomyopathy. Patients were excluded if cardiogenic shock was secondary to cardiotoxic poisoning, if they were pregnant, or were under legal protective supervision (e.g., court-appointed guardianship).

Data Collection

Clinical, laboratory and echocardiographic parameters were prospectively collected at three predefined timepoints:Day 0 (ECMO implantation), Day 3–5 (early course) and Day of ECMO weaning (ECMO explantation). Blood samples were drawn at each timepoint, immediately centrifuged and plasma was stored at –80°C until analysis. Samples were then transferred to the XXX for cells analysis and INSERM UMR942 laboratory (Paris, France) for cytokine analysis. Quantification of adrenoreceptors β1 and β2 on monocytes and CD4⁺ T lymphocytes was performed using flow cytometry. Cytokine concentrations in plasma, including interleukin (IL) -33, IL-1β, IL-2, IL-4, IL-5, IL-6, IL-8/CXCL8, IL-10, GDF-15, IL-12p70, IL-17A/CTLA-8, IFN-γ, TNF-α and granzyme B, were measured using [XXX].

Assessment of ADRB1 and ADRB2

The expression of β1 (ADRB1) and β2 (ADRB2) adrenergic receptors on circulating immune cell, namely CD4+ T lymphocytes and monocytes, was assessed via flow cytometry. XXX [to be completed from Manon]. Results were expressed as the percentage of cells that were positive for the respective expression receptors.

##Cytokine Quantification

The aforementioned panel of circulating cytokines and immune mediators was quantified from plasma samples. These biomarkers were selected for their roles in pro- and anti-inflammatory signaling. Measurements were performed using a multiplex bead-based immunoassay (Luminex platform), following the manufacturer’s instructions. XXX [to be completed from Malha].

##ECMO Weaning Protocol

A standardized weaning protocol was applied to all patients prior to ECMO explantation. At a constant mean arterial pressure (MAP) between 65 and 75 mmHg, cardiac output was assessed by transthoracic echocardiography at two predefined timepoints: 1. Baseline under full support (ECLS flow at 3.0 L/min) and 2. Low-flow test after 45 minutes of reduced ECMO flow (1.5 L/min). Measurements included velocity time integral (VTI) at the left ventricular outflow tract and left ventricular ejection fraction. Changes in norepinephrine and dobutamin dosage during this period were also documented. The decision to explant ECMO was based on clinical tolerance, hemodynamic stability and echocardiographic parameters during the weaning trial.

#Data load

General consideration

# Blood sample at day 0 before ECMO implant?

df %>%
  select(J0_date, j0_bilan_date) %>%
  distinct() %>%
  slice(1:10)

# A tibble: 10 × 2
   J0_date    j0_bilan_date
   <chr>      <chr>        
 1 16.10.2017 17.10.2017   
 2 25.10.2017 25.10.2017   
 3 17.01.2018 17.01.2018   
 4 25.01.2018 25.01.2018   
 5 06.02.2018 06.02.2018   
 6 20.02.2018 20.02.2018   
 7 02.03.2018 02.03.2018   
 8 30.05.2018 30.05.2018   
 9 14.06.2018 14.06.2018   
10 27.06.2018 27.06.2018   

df <- df %>%
  mutate(
    J0_date = as.Date(J0_date, format = "%d.%m.%Y"),
    j0_bilan_date = as.Date(j0_bilan_date, format = "%d.%m.%Y"),
    bilan_after_ecmo = case_when(
      is.na(J0_date) | is.na(j0_bilan_date) ~ NA_character_,
      j0_bilan_date > J0_date ~ "Yes",
      TRUE ~ "No"
    )
  )

table(df$bilan_after_ecmo, useNA = "ifany")

 No Yes 
 48   1 

## 2 patients had blood sample AFTER ECMO implant.

df %>%
  filter(bilan_after_ecmo == "Yes") %>%
  select(ID, j0_bilan_date, J0_date, bilan_after_ecmo)

# A tibble: 1 × 4
  ID    j0_bilan_date J0_date    bilan_after_ecmo
  <chr> <date>        <date>     <chr>           
1 1-FG  2017-10-17    2017-10-16 Yes             

# Patient 1 (2017-10-17 > 2017-10-16) and 15 (2019-04-19 > 2018-04-19). Patients 15 has probably an error. In PAtinet 1 the blood sample was after 12 hours.


df <- df %>%
  # Applica solo alle righe dove bilan_after_ecmo è "No"
  mutate(
    heure_check = case_when(
      bilan_after_ecmo == "No" &
        !is.na(J0_heure) & !is.na(j0_bilan_heure) ~ if_else(
          j0_bilan_heure > J0_heure, "Bilan after ECMO (same day)",
          "Bilan before/same as ECMO (same day)"
        ),
      TRUE ~ NA_character_
    )
  )

table(df$bilan_after_ecmo, useNA = "ifany")

 No Yes 
 48   1 

df %>%
  filter(heure_check == "Bilan after ECMO (same day)") %>%
  select(ID, J0_date, J0_heure, j0_bilan_date, j0_bilan_heure, heure_check)

# A tibble: 0 × 6
# ℹ 6 variables: ID <chr>, J0_date <date>, J0_heure <time>,
#   j0_bilan_date <date>, j0_bilan_heure <time>, heure_check <chr>

# Among patients with laboratory tests on the same day of ECMO initiation, none had blood sampling performed after ECMO implantation.



# Mean time between blood sample and ECMO implant, patients 1 e 15 excluded.

df_clean <- df %>%
  filter(!str_starts(ID, "1-") & !str_starts(ID, "15-")) %>%
  mutate(
    j0_bilan_date = as.Date(j0_bilan_date, format = "%d.%m.%Y"),
    J0_date = as.Date(J0_date, format = "%d.%m.%Y"),
    datetime_bilan = as.POSIXct(paste(j0_bilan_date, as.character(j0_bilan_heure)), format = "%Y-%m-%d %H:%M", tz = "UTC"),
    datetime_ecmo  = as.POSIXct(paste(J0_date, as.character(J0_heure)), format = "%Y-%m-%d %H:%M", tz = "UTC"),
    time_diff_mins = as.numeric(difftime(datetime_bilan, datetime_ecmo, units = "mins"))
  ) %>%
  filter(!is.na(time_diff_mins))

df_clean %>%
  summarise(
    mean_minutes = mean(time_diff_mins),
    sd_minutes = sd(time_diff_mins),
    min = min(time_diff_mins),
    max = max(time_diff_mins)
  )

# A tibble: 1 × 4
  mean_minutes sd_minutes   min   max
         <dbl>      <dbl> <dbl> <dbl>
1        -205.       145.  -660     0

In 2 patients, the day 0 laboratory test occurred after ECMO initiation. Day-0 blood samples were collected on average 3.4 hours (SD 2.4 hours) before ECMO implantation, ranging from 11 hours before to the exact time of ECMO cannulation.

DATA Management

df <- df %>%
  mutate(
    # Basic formatting
    Temp_D0 = round(Temp_D0, 1),
    Age = as.numeric(Age),
    Gender = case_when(
      Gender == 0 ~ "Men",
      Gender == 1 ~ "Women",
      TRUE ~ as.character(Gender)
    ),
    
    Preecmo_tte_diam_LV = Preecmo_tte_diam_LV * 10,
    
    VIS_D0 = (DOBUcum_D0 * 1000) / (Weight * 1440) + 
         100 * (NADcum_D0 * 1000) / (Weight * 1440),

    # Outcome recoding
    Outcome = case_when(
      Outcome == 1 ~ "ECMO Weaning",
      Outcome == 2 ~ "Bridge to Transplant",
      Outcome == 3 ~ "Bridge to LVAD",
      Outcome == 4 ~ "Death",
      TRUE ~ as.character(Outcome)
    ),

    Pulsepressure_D0 = PAS_D0 - PAD_D0,
    
    Outcome_censored = as.numeric(case_when(
      Outcome == "Death" ~ "1",
      Outcome %in% c("ECMO Weaning", "Bridge to Transplant", "Bridge to LVAD") ~ "0"
    )),

    outcome_date = str_replace_all(outcome_date, "\\.", "/"),
    icu_admiss_date = str_replace_all(icu_admiss_date, "\\.", "/"),
    icu_discharge_date = str_replace_all(icu_discharge_date, "\\.", "/"),

    outcome_date = dmy(outcome_date),
    icu_admiss_date = dmy(icu_admiss_date),
    icu_discharge_date = dmy(icu_discharge_date),
    date_levo = dmy(date_levo),
    ecmo_start_date = dmy(ecmo_start_date),
    ecmo_stop_date = dmy(ecmo_stop_date),

    diff_days = time_length(interval(icu_admiss_date, outcome_date), "days"),
    Time_hosp = time_length(interval(icu_admiss_date, icu_discharge_date), "days"),
    as.numeric(icu_discharge_date - icu_admiss_date),
    Time_ecmo = time_length(interval(ecmo_start_date, ecmo_stop_date), "days"),
    Time_levo = time_length(interval(ecmo_start_date, date_levo), "days"),

    diff_days_J28 = case_when(
      diff_days > 28 ~ 28,
      TRUE ~ diff_days
    ),

    Outcome_censored_28 = case_when(
      diff_days > 28 & Outcome_censored == 1 ~ 0,
      TRUE ~ Outcome_censored
    ),

    ADR1_D0_median = case_when(
      ADR1_D0 > median(ADR1_D0, na.rm = TRUE) ~ "High",
      ADR1_D0 <= median(ADR1_D0, na.rm = TRUE) ~ "Low"
    ),

    ADR2_D0_median = case_when(
      ADR2_D0 > median(ADR2_D0, na.rm = TRUE) ~ "High",
      ADR2_D0 <= median(ADR2_D0, na.rm = TRUE) ~ "Low"
    ),
    
    ADR1m_D0_median = case_when(
      ADR1m_J0 > median(ADR1m_J0, na.rm = TRUE) ~ "High",
      ADR1m_J0 <= median(ADR1m_J0, na.rm = TRUE) ~ "Low",
      TRUE ~ NA_character_
    ),
    ADR2m_D0_median = case_when(
      ADR2m_J0 > median(ADR2m_J0, na.rm = TRUE) ~ "High",
      ADR2m_J0 <= median(ADR2m_J0, na.rm = TRUE) ~ "Low",
      TRUE ~ NA_character_
    ),

    Betablocker = case_when(
      Betablocker == 1 ~ "YES",
      Betablocker == 0 ~ "NO"
    ),

    dose_tot_dobu_yes_no = droplevels(as.factor(case_when(
      dose_tot_dobu > 0 ~ "YES",
      dose_tot_dobu == 0 ~ "NO"
    ))),

    levo_in_icu = case_when(
      Levosimendan == 1 & Time_levo >= 0 ~ "Yes",  
      Levosimendan == 1 & Time_levo < 0 ~ "No",
      Levosimendan == 0 ~ "No"
    ),

    Levosimendan = case_when(
      Levosimendan == 1 ~ "YES",
      Levosimendan == 0 ~ "NO",
      TRUE ~ NA_character_ 
    )
  )

var.cat <- c("Gender", "HTN", "DM", "Immunodepression", "CKD", "Neurologic_deficit", "Chronic_respiratory_disease", "HF", "Cardiac_arrest_before_canul", "IABP_D0", "Intubation", "EER_D0", "Outcome", "Betablocker", "Alphablocker", "Levosimendan", "Death")

var.cont <- c("Age", "BMI", "PAS_D0", "PAD_D0", "PAM_D0", "Pulsepressure_D0", "HR_D0", "Temp_D0", "Preecmo_tte_ef", "Preecmo_tte_vtiao", "Preecmo_tte_diam_LV", "Preecmo_tte_tapse", "Preecmo_tte_vmax_s_mit_lat_LV", "NADcum_D0", "DOBUcum_D0", "VIS_D0", "pH_D0", "Sao2_D0", "Paco2_D0", "Lact_D0", "Hco3_D0", "Creat_D0", "Urea_D0", "ALAT_D0", "ASAT_D0", "Bili_D0", "PT_D0", "Tropo_i_hs_D0", "Ntprobnp_D0", "DPP3_D0", "Plq_D0", "Lenght_eer", "Lenght_vm", "Time_hosp", "Time_ecmo")

var.tot <- c(var.cat, var.cont)

TABLE 1: Characteristics of the population

table1

Table 1: Characteristics of the population

Variables	N	Overall
Age (years)	50	59.00 [49.00, 66.00]
Gender = Female (%)	50	8 (16.3)
Body Mass Index (kg/m²)	47	26.98 [23.72, 30.00]
Hypertension (%)	50	19 (38.8)
Diabetes Mellitus (%)	50	10 (20.4)
Heart Failure (%)	50	17 (34.7)
Chronic Kidney Disease (%)	50	5 (10.2)
Chronic Respiratory Disease (%)	50	3 ( 6.1)
Neurological Deficit (%)	50	4 ( 8.2)
Immunosuppression (%)	50	4 ( 8.2)
Systolic Blood Pressure (mmHg)	50	96.00 [86.00, 109.00]
Diastolic Blood Pressure (mmHg)	50	63.00 [56.00, 72.00]
Mean Arterial Pressure (mmHg)	50	72.00 [68.00, 81.00]
Pulse Pressure (mmHg)	50	31.00 [22.00, 47.00]
Heart Rate (bpm)	50	96.00 [86.00, 118.00]
Temperature (°C)	50	36.00 [35.30, 36.80]
Cardiac Arrest Before Cannulation (%)	50	27 (55.1)
Intra-Aortic Balloon Pump (%)	50	23 (46.9)
Intubation (%)	50	36 (73.5)
Renal Replacement Therapy (%)	50	7 (14.3)
Pre-ECMO Ejection Fraction (%)	46	20.00 [15.00, 25.00]
Pre-ECMO Aortic Velocity Time Integral (cm)	41	7.00 [5.27, 8.25]
Pre-ECMO Left Ventricular Diameter (mm)	31	59.00 [48.50, 65.25]
Pre-ECMO TAPSE (mm)	38	12.50 [8.00, 16.00]
Pre-ECMO Mitral Lateral Velocity (cm/s)	26	4.50 [4.00, 5.80]
Beta-blocker (%)	50	22 (44.9)
Alpha-blocker (%)	50	3 ( 6.1)
Cumulative Norepinephrine Dose (mg)	50	20.00 [0.00, 86.50]
Cumulative Dobutamine Dose (mg)	50	200.00 [0.00, 510.00]
Vasoactive-Inotropic Score (VIS) at Day 0	48	18.52 [4.86, 58.75]
Levosimendan (%)	50	29 (59.2)
pH	50	7.42 [7.31, 7.49]
Oxygen Saturation (%)	49	98.00 [96.00, 99.00]
PaCO2 (mmHg)	50	32.70 [25.90, 37.60]
Lactate (mmol/L)	49	2.75 [1.70, 4.70]
Bicarbonate (mmol/L)	50	19.00 [16.80, 23.00]
Creatinine (µmol/L)	50	149.00 [104.00, 202.00]
Urea (mg/dL)	49	11.21 [8.14, 18.86]
ALT (IU/L)	50	264.00 [91.00, 804.00]
AST (IU/L)	50	656.00 [171.00, 1202.00]
Bilirubin (mg/dL)	49	21.00 [12.00, 38.50]
Prothrombin Time (s)	38	56.00 [37.00, 73.00]
High Sensitivity Troponin I (ng/L)	35	34363.50 [383.25, 186316.00]
NT-proBNP (pg/mL)	24	10882.00 [2248.00, 26995.00]
DPP3 (ng/mL)	50	72.90 [50.00, 164.10]
Platelet Count (10⁹/L)	50	176.00 [118.00, 271.00]
Outcome (%)	50
Death in ICU	50	22 (44.9)
Bridge to LVAD	50	4 ( 8.2)
Bridge to Transplant	50	2 ( 4.1)
ECMO Weaning	50	21 (42.9)
Duration of Renal Replacement Therapy (days)	50	0.00 [0.00, 5.00]
Duration of Mechanical Ventilation (days)	50	5.00 [2.00, 10.00]
ECMO Duration (days)	50	6.00 [3.00, 9.00]
ICU Length of Stay (days)	50	14.00 [6.00, 19.00]
28-day Mortality (%)	50	28 (57.1)
	50	49

Statistical analysis

Baseline clinical, demographic and hemodynamic characteristics were summarized for the overall study population (n = 50). One patient was excluded from the analysis due to being under legal guardianship. Categorical variables were reported as counts and percentages, while continuous variables were expressed as medians and interquartile ranges (IQR). The normality of continuous variables was assessed visually using histograms.

Check for normality

for (var in var.cont) {
  if (var %in% colnames(T1)) {
    print(
      ggplot(T1, aes_string(var)) +
        geom_histogram(bins = 15, fill = "blue", alpha = 0.5) +
        ggtitle(paste("Distribuzione di", var))
    )
  } else {
    message(paste("❌ Variabile", var, "non trovata in T1"))
  }
}

Table 2a with hemodynamic parameters at weaning day

table2a

# A tibble: 8 × 2
  Parameter                     `Median [IQR]`     
  <chr>                         <chr>              
1 ECMO flow (l)                 2.4 [1.6–3], n=33  
2 Dobutamine cumulative (mg)    0 [0–50], n=33     
3 Noradrenaline cumulative (mg) 0 [0–6.2], n=33    
4 MAP_baseline (mmHg)           75 [72–80], n=33   
5 HR (bpm)                      100 [84–103], n=33 
6 FEVG (%)                      35 [30–43.2], n=32 
7 Mitral S' (cm/s)              7.4 [6–9], n=31    
8 Lactate (mmol/L)              1.1 [0.9–1.4], n=32

Table 2b hemodynamic parameters during the weaning test

# Create new variables for flow condition and unify variable names

df_long <- df %>%
  select(
    ID,
    jsev_epr_sev_t0_3l_itv,
    jsev_epr_sev_t0_3l_pam,
    jsev_epr_sev_t0_3l_NAd,
    `jsev_epr_sev_t1_1,5l_itv`,
    `jsev_epr_sev_t1_1,5l_pam`,
    `jsev_epr_sev_t1_1,5l_NAd`,
    Weight,
  ) %>%
  mutate(jsev_epr_sev_t0_3l_NAd_w = jsev_epr_sev_t0_3l_NAd/Weight,
         `jsev_epr_sev_t1_1,5l_NAd_w` = `jsev_epr_sev_t1_1,5l_NAd`/Weight)%>%
  rename(
    `VTI (cm)` = jsev_epr_sev_t0_3l_itv,
    `MAP (mmHg)` = jsev_epr_sev_t0_3l_pam,
    `Noradrenalin (µg/kg/h)` = jsev_epr_sev_t0_3l_NAd_w,
    `VTI_1.5L (cm)` = `jsev_epr_sev_t1_1,5l_itv`,
    `MAP_1.5L (mmHg)` = `jsev_epr_sev_t1_1,5l_pam`,
    `Noradrenalin_1.5L (µg/kg/h)` = `jsev_epr_sev_t1_1,5l_NAd_w`
  ) %>%
  pivot_longer(
    cols = -c(ID,jsev_epr_sev_t0_3l_NAd,`jsev_epr_sev_t1_1,5l_NAd`,Weight),
    names_to = "Parameter_raw",
    values_to = "Value"
  ) %>%
  mutate(
    Value = as.numeric(Value),
    Condition = case_when(
      str_detect(Parameter_raw, "1.5L") ~ "1.5L",
      TRUE ~ "3L"
    ),
    Condition = factor(Condition, levels = c("3L", "1.5L")),
    Parameter = str_replace(Parameter_raw, "_1.5L", "")
  )


# Create summary table

table2b <- df_long %>%
  group_by(Parameter, Condition) %>%
  summarise(
    n = sum(!is.na(Value)),
    Median = round(median(Value, na.rm = TRUE), 2),
    Q1 = round(quantile(Value, 0.25, na.rm = TRUE), 2),
    Q3 = round(quantile(Value, 0.75, na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  mutate(
    Value_IQR = paste0(Median, " [", Q1, "–", Q3, "]; n=", n)
  ) %>%
  select(Parameter, Condition, Value_IQR) %>%
  pivot_wider(names_from = Condition, values_from = Value_IQR)

table2b

# A tibble: 3 × 3
  Parameter              `3L`                    `1.5L`              
  <chr>                  <chr>                   <chr>               
1 MAP (mmHg)             75.5 [70.75–83.5]; n=20 72.5 [70–80.5]; n=28
2 Noradrenalin (µg/kg/h) 0 [0–0]; n=22           0 [0–0]; n=29       
3 VTI (cm)               13 [12–15]; n=21        14 [12–16]; n=29    

# Wilcoxon test for a parameter

wilcox.test(Value ~ Condition,
            data = df_long %>% filter(Parameter == "VTI (cm)"))

    Wilcoxon rank sum test with continuity correction

data:  Value by Condition
W = 271.5, p-value = 0.5204
alternative hypothesis: true location shift is not equal to 0

No significant difference between VTI at 3L and 1.5L flow conditions at weaning day (Wilcoxon test, p = 0.5).

Figure 1: Association between ADRB expressions and outcomes

F1 <- df %>%
  select(ID,ADR1_J0, ADR1_J3_J5, ADR1_JS, Outcome)%>%
  mutate(Outcomes = factor(case_when(Outcome=="Bridge to LVAD"|Outcome=="Bridge to Transplant"~ "Bridge to transplant or LVAD",T~Outcome), levels=c("Death","Bridge to transplant or LVAD","ECMO Weaning")))%>%
  pivot_longer(
    cols = c(ADR1_J0, ADR1_J3_J5, ADR1_JS), 
    names_to = "ADR1_time", 
    values_to = "value"
  )%>%
  mutate(time = as.numeric(case_when(ADR1_time=="ADR1_J0"~1,
                          ADR1_time=="ADR1_J3_J5"~2,
                          ADR1_time=="ADR1_JS"~3)),
         value_log = log10(value))

count_data <- F1 %>%
  filter(!is.na(value))%>%
  group_by(ADR1_time, Outcomes) %>%
  summarise(n = n(), .groups = "drop")%>%
  filter(!is.na(Outcomes))


# ===========================
# Modélisation
# ===========================

# Modèle avec interaction
mod1 <- lmerTest::lmer(
  value_log ~ time * Outcomes + (1 + time | ID),
  data = F1
)


# Résumé du modèle
summary(mod1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: value_log ~ time * Outcomes + (1 + time | ID)
   Data: F1

REML criterion at convergence: 100.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.35927 -0.53586  0.03986  0.52500  1.90217 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 ID       (Intercept) 0.11537  0.3397        
          time        0.01505  0.1227   -0.71
 Residual             0.09453  0.3075        
Number of obs: 93, groups:  ID, 46

Fixed effects:
                                          Estimate Std. Error       df t value
(Intercept)                               -0.15452    0.16111 61.94485  -0.959
time                                       0.12703    0.08955 54.66302   1.419
OutcomesBridge to transplant or LVAD       0.09504    0.28976 38.98471   0.328
OutcomesECMO Weaning                       0.32514    0.21202 48.22743   1.534
time:OutcomesBridge to transplant or LVAD  0.09162    0.13978 34.93053   0.655
time:OutcomesECMO Weaning                 -0.07095    0.10946 43.99886  -0.648
                                          Pr(>|t|)
(Intercept)                                  0.341
time                                         0.162
OutcomesBridge to transplant or LVAD         0.745
OutcomesECMO Weaning                         0.132
time:OutcomesBridge to transplant or LVAD    0.516
time:OutcomesECMO Weaning                    0.520

Correlation of Fixed Effects:
            (Intr) time   OBttoL OECMOW t:ttoL
time        -0.877                            
OtcBttoLVAD -0.556  0.487                     
OtcmsECMOWn -0.760  0.666  0.422              
t:OBttoLVAD  0.562 -0.641 -0.860 -0.427       
tm:OtcECMOW  0.717 -0.818 -0.399 -0.860  0.524

# Modèle sans interaction (modèle réduit)
mod_no_interaction <- lmerTest::lmer(
  value_log ~ time + Outcomes + (1 + time | ID),
  data = F1
)

# Comparaison des modèles (test de l'interaction)
interact <- anova(mod_no_interaction, mod1)
# Interprétation : p = 0.58 → interaction non significative

# ===========================
# Tests des effets fixes
# ===========================

# ANOVA avec approximation de Kenward-Roger pour les p-values
a <- anova(mod_no_interaction, ddf = "Kenward-Roger")
print(a)

Type III Analysis of Variance Table with Kenward-Roger's method
          Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
time     0.48396 0.48396     1 32.800  5.0378 0.03167 *
Outcomes 0.46031 0.23015     2 36.375  2.3938 0.10552  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ===========================
# Vérification des hypothèses
# ===========================

# 1. Diagrammes de diagnostic généraux
plot(mod1)  # résidus vs prédictions par groupe

# 2. Résidus vs valeurs ajustées (homoscédasticité)
plot(fitted(mod1), residuals(mod1),
     xlab = "Valeurs ajustées", ylab = "Résidus",
     main = "Résidus vs ajustés")

# 3. Normalité des résidus
qqnorm(residuals(mod1))
qqline(residuals(mod1), col = "blue")

# 4. Autocorrélation des résidus (si données longitudinales)
acf(residuals(mod1), main = "ACF des résidus")

# 5. Normalité des pentes aléatoires (effets aléatoires sur time)
qqnorm(ranef(mod1)$ID[, "time"],
       main = "QQ-plot des pentes aléatoires (time)")
qqline(ranef(mod1)$ID[, "time"], col = "blue")

#extrait P values
annot_pvals <- data.frame(
  ADR1_time = "ADR1_J0",
  value_log = max(F1$value_log, na.rm = TRUE) * 1.1,
  label = "p[time] == 0.024 * '\n' * p[outcomes] == 0.066"
)

# ===========================
# Figure
# ===========================
Figure_1 <- ggplot(F1, aes(x = ADR1_time, y = value_log, fill = Outcomes)) +
  geom_boxplot(position = position_dodge(0.7), width = 0.6) +  
  scale_y_log10() +  
  labs(
    title = "Association between ADRB1 and outcomes",
    x = "times of measurement",
    y = "Percentage of T4 lymphocyte expressing ADRB1\n (log 10 scale)",
  ) +
  geom_text(data = count_data, aes(x = ADR1_time, y =- min(F1$value_log, na.rm = TRUE) *0.1, 
                                   label = paste0("n=", n), group = Outcomes), 
            position = position_dodge(0.7), size = 4, vjust = 1)+
    scale_x_discrete(labels = c("ADR1_J0" = "implantation", "ADR1_J3_J5" = "day 3 to 5", "ADR1_JS" = "explantation")) +
  scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",                   
    "Bridge to transplant or LVAD" = "#377EB8",   
    "Death" = "#F8766D"                           
  )) +
  theme(legend.position = "none")+ 
  annotate(
    "text",
    x = 1,  # position sur l’axe x : 1 = "ADR1_J0"
    y = max(F1$value_log, na.rm = TRUE) * 1.1,  # position en haut du graphique
    label = "p[time] == 0.024*';'~p[outcomes] == 0.066",
    parse = TRUE,
    size = 4,
    hjust = 0
  )

#####################
## ADRB2 
#####################


F1_1 <- df %>%
  select(ID,ADR2_J0, ADR2_J3_J5, ADR2_JS, Outcome)%>%
  mutate(Outcomes = factor(case_when(Outcome=="Bridge to LVAD"|Outcome=="Bridge to Transplant"~ "Bridge to transplant or LVAD",T~Outcome), levels=c("Death","Bridge to transplant or LVAD","ECMO Weaning")))%>%
  pivot_longer(
    cols = c(ADR2_J0, ADR2_J3_J5, ADR2_JS), 
    names_to = "ADR2_time", 
    values_to = "value"
  )%>%
  mutate(time = as.numeric(case_when(ADR2_time=="ADR2_J0"~1,
                          ADR2_time=="ADR2_J3_J5"~2,
                          ADR2_time=="ADR2_JS"~3)),
         value_log = log10(value))


count_data <- F1_1 %>%
  filter(!is.na(value))%>%
  group_by(ADR2_time, Outcomes) %>%
  summarise(n = n(), .groups = "drop")%>%
  filter(!is.na(Outcomes))

# ===========================
# Modélisation
# ===========================

# Modèle avec interaction
mod1 <- lmerTest::lmer(
  value_log ~ time * Outcomes + (1 + time | ID),
  data = F1_1
)

# Résumé du modèle
summary(mod1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: value_log ~ time * Outcomes + (1 + time | ID)
   Data: F1_1

REML criterion at convergence: 101.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.25460 -0.42375 -0.02913  0.60839  1.65730 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 ID       (Intercept) 0.15994  0.3999        
          time        0.02624  0.1620   -0.77
 Residual             0.08183  0.2861        
Number of obs: 95, groups:  ID, 46

Fixed effects:
                                          Estimate Std. Error       df t value
(Intercept)                               -0.02846    0.15731 49.17843  -0.181
time                                       0.08396    0.08554 43.13429   0.982
OutcomesBridge to transplant or LVAD       0.08591    0.29201 34.60684   0.294
OutcomesECMO Weaning                       0.26694    0.20966 41.35055   1.273
time:OutcomesBridge to transplant or LVAD  0.12482    0.14021 28.85470   0.890
time:OutcomesECMO Weaning                 -0.04361    0.10803 35.41604  -0.404
                                          Pr(>|t|)
(Intercept)                                  0.857
time                                         0.332
OutcomesBridge to transplant or LVAD         0.770
OutcomesECMO Weaning                         0.210
time:OutcomesBridge to transplant or LVAD    0.381
time:OutcomesECMO Weaning                    0.689

Correlation of Fixed Effects:
            (Intr) time   OBttoL OECMOW t:ttoL
time        -0.868                            
OtcBttoLVAD -0.539  0.468                     
OtcmsECMOWn -0.750  0.651  0.404              
t:OBttoLVAD  0.530 -0.610 -0.860 -0.397       
tm:OtcECMOW  0.688 -0.792 -0.370 -0.858  0.483

# Modèle sans interaction (modèle réduit)
mod_no_interaction <- lmerTest::lmer(
  value_log ~ time + Outcomes + (1 + time | ID),
  data = F1_1
)

# Comparaison des modèles (test de l'interaction)
interact <- anova(mod_no_interaction, mod1)
# Interprétation : p = 0.76 → interaction non significative

# ===========================
# Tests des effets fixes
# ===========================

# ANOVA avec approximation de Kenward-Roger pour les p-values
a <- anova(mod_no_interaction, ddf = "Kenward-Roger")
print(a)

Type III Analysis of Variance Table with Kenward-Roger's method
          Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
time     0.25391 0.25391     1 33.452  3.0552 0.08965 .
Outcomes 0.45728 0.22864     2 36.789  2.7496 0.07711 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ===========================
# Vérification des hypothèses
# ===========================

# 1. Diagrammes de diagnostic généraux
plot(mod1)  # résidus vs prédictions par groupe

# 2. Résidus vs valeurs ajustées (homoscédasticité)
plot(fitted(mod1), residuals(mod1),
     xlab = "Valeurs ajustées", ylab = "Résidus",
     main = "Résidus vs ajustés")

# 3. Normalité des résidus
qqnorm(residuals(mod1))
qqline(residuals(mod1), col = "blue")

# 4. Autocorrélation des résidus (si données longitudinales)
acf(residuals(mod1), main = "ACF des résidus")

# 5. Normalité des pentes aléatoires (effets aléatoires sur time)
qqnorm(ranef(mod1)$ID[, "time"],
       main = "QQ-plot des pentes aléatoires (time)")
qqline(ranef(mod1)$ID[, "time"], col = "blue")

#extrait P values
annot_pvals <- data.frame(
  ADR2_time = "ADR2_J0",
  value_log = max(F1_1$value_log, na.rm = TRUE) * 1.1,
  label = "p[time] == 0.07 * '\n' * p[outcomes] == 0.05"
)

Figure_1.1 <- ggplot(F1_1, aes(x = ADR2_time, y = value_log, fill = Outcomes)) +
  geom_boxplot(position = position_dodge(0.7), width = 0.6) +  
  scale_y_log10() +  
  labs(
    title = "Association between ADRB2 and outcomes",
    x = "times of measurement",
    y = "Percentage of T4 lymphocyte expressing ADRB2\n (log 10 scale)",
  ) +
  geom_text(data = count_data, aes(x = ADR2_time, y = -min(F1_1$value_log, na.rm = TRUE) * 0.01, 
                                   label = paste0("n=", n), group = Outcomes), 
            position = position_dodge(0.7), size = 4, vjust = 1)+
      scale_x_discrete(labels = c("ADR2_J0" = "implantation", "ADR2_J3_J5" = "day 3 to 5", "ADR2_JS" = "explantation"))+
    scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",                   
    "Bridge to transplant or LVAD" = "#377EB8",   
    "Death" = "#F8766D"                           
  )) +
    annotate(
    "text",
    x = 1,  # position sur l’axe x : 1 = "ADR1_J0"
    y = max(F1_1$value_log, na.rm = TRUE) * 1.1,  # position en haut du graphique
    label = "p[time] == 0.07*';'~p[outcomes] == 0.05",
    parse = TRUE,
    size = 4,
    hjust = 0
  )



Fig_1 <- Figure_1 #+ theme(plot.title = element_blank())
Fig_1.1 <- Figure_1.1 #+ theme(plot.title = element_blank())

Figure1_tot <- Fig_1 + Fig_1.1 +
  plot_annotation(title = "Association between percentage of T4 lymphocyte expressing ADRB and Outcome")

Figure1_tot

Figure 1m: Association between ADRB (monocytes) and outcomes

# Prepare data
F1m <- df %>%
  select(ID, ADR1m_J0, ADR1m_J3_J5, ADR1m_JS, Outcome) %>%
  mutate(Outcomes = factor(case_when(
    Outcome == "Bridge to LVAD" | Outcome == "Bridge to Transplant" ~ "Bridge to transplant or LVAD",
    TRUE ~ Outcome
  ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))) %>%
  pivot_longer(
    cols = c(ADR1m_J0, ADR1m_J3_J5, ADR1m_JS), 
    names_to = "ADR1_time", 
    values_to = "value"
  ) %>%
  mutate(
    time = case_when(
      ADR1_time == "ADR1m_J0" ~ 1,
      ADR1_time == "ADR1m_J3_J5" ~ 2,
      ADR1_time == "ADR1m_JS" ~ 3
    )
  ) %>%
  filter(!is.na(value), is.finite(value))

F1m <- F1m %>%
  mutate(value_log = log10(value)) %>%
  filter(!is.na(value_log), is.finite(value_log))  


# Count data for annotation
count_data_m <- F1m %>%
  filter(!is.na(value)) %>%
  group_by(ADR1_time, Outcomes) %>%
  summarise(n = n(), .groups = "drop") %>%
  filter(!is.na(Outcomes))

# ===========================
# Mixed effects model
# ===========================
# Model with interaction
mod1_m <- lmerTest::lmer(
  value_log ~ time * Outcomes + (1 + time | ID),
  data = F1m
)

# Model without interaction
mod_no_interaction_m <- lmerTest::lmer(
  value_log ~ time + Outcomes + (1 + time | ID),
  data = F1m
)

# Interaction test
interact_m <- anova(mod_no_interaction_m, mod1_m)
# Interpretation : p = 0.64 → interaction not significant

# ANOVA for fixed effect (withouth interaction)
a_m <- anova(mod_no_interaction_m, ddf = "Kenward-Roger")
print(a_m)

Type III Analysis of Variance Table with Kenward-Roger's method
           Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
time     0.056058 0.056058     1 33.493  0.3582 0.5535
Outcomes 0.105353 0.052676     2 39.504  0.3365 0.7163

# ===========================
# Diagnostic plots
# ===========================
# 1. Residual diagnostic plots
plot(mod1_m)

# 2. Residuals vs Fitted
plot(fitted(mod1_m), residuals(mod1_m),
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs Fitted")

# 3. Normality of residuals
qqnorm(residuals(mod1_m))
qqline(residuals(mod1_m), col = "blue")

# 4. Autocorrelation of residuals
acf(residuals(mod1_m), main = "ACF of residuals")

# 5. Normality of random intercepts (optional)
qqnorm(ranef(mod1_m)$ID[, "(Intercept)"],
       main = "QQ plot of random intercepts (ID)")
qqline(ranef(mod1_m)$ID[, "(Intercept)"], col = "blue")

# ===========================
# Figure
# ===========================
Figure_1m <- ggplot(F1m, aes(
    x = ADR1_time,
    y = value,
    fill = Outcomes,
    group = interaction(ADR1_time, Outcomes)  
)) +
  geom_boxplot(position = position_dodge(0.7), width = 0.6) +
  scale_y_log10() +
  labs(
    title = "Association between ADRB1 expression in monocytes and outcomes",
    x = "Time of measurement",
    y = "Percentage of monocytes expressing ADRB1\n(log10 scale)"
  ) +
  geom_text(data = count_data_m, aes(
    x = ADR1_time,
    y = min(F1m$value[F1m$value > 0], na.rm = TRUE) * 0.5,
    label = paste0("n=", n),
    group = interaction(ADR1_time, Outcomes)
  ),
  position = position_dodge(0.7),
  size = 4,
  vjust = 1) +
  scale_x_discrete(labels = c(
    "ADR1m_J0" = "implantation",
    "ADR1m_J3_J5" = "day 3 to 5",
    "ADR1m_JS" = "explantation"
  )) +
    scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",                   
    "Bridge to transplant or LVAD" = "#377EB8",   
    "Death" = "#F8766D"                           
  )) +
  theme(legend.position = "none") +
  annotate(
    "text",
    x = 1,
    y = max(F1m$value, na.rm = TRUE) * 1.1,
    label = "p[time] == 0.955*';'~p[outcomes] == 0.649",
    parse = TRUE,
    size = 4,
    hjust = 0
  )

#####################
## ADRB2 (monocytes)
#####################

# ==========================
# Prepare data (non log)
# ==========================
F1m_2 <- df %>%
  select(ID, ADR2m_J0, ADR2m_J3_J5, ADR2m_JS, Outcome) %>%
  mutate(
    Outcomes = factor(case_when(
      Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
      TRUE ~ Outcome
    ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
  ) %>%
  pivot_longer(
    cols = c(ADR2m_J0, ADR2m_J3_J5, ADR2m_JS),
    names_to = "ADR2_time",
    values_to = "value"
  ) %>%
  mutate(
    time = case_when(
      ADR2_time == "ADR2m_J0" ~ 1,
      ADR2_time == "ADR2m_J3_J5" ~ 2,
      ADR2_time == "ADR2m_JS" ~ 3
    )
  ) %>%
  filter(!is.na(value), is.finite(value))

F1m_2 <- F1m_2 %>%
  mutate(value_log = log10(value)) %>%
  filter(!is.na(value_log), is.finite(value_log))

# ==========================
# Count available data
# ==========================
count_data_m2 <- F1m_2 %>%
  group_by(ADR2_time, Outcomes) %>%
  summarise(n = n(), .groups = "drop")

# ==========================
# Mixed effects model (simplified)
# ==========================
mod1_m2 <- lmerTest::lmer(
  value_log ~ time * Outcomes + (1 + time | ID),
  data = F1m_2
)


mod_no_interaction_m2 <- lmerTest::lmer(
  value_log ~ time + Outcomes + (1 + time | ID),
  data = F1m_2
)

# Test interaction
interact_m2 <- anova(mod_no_interaction_m2, mod1_m2)
# Interpretation : p = 0.76 → interaction not significant

# Fixed effect ANOVA (no interaction)
a_m2 <- anova(mod_no_interaction_m2, ddf = "Kenward-Roger")
print(a_m2)

Type III Analysis of Variance Table with Kenward-Roger's method
          Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
time     0.33421 0.33421     1 32.054  1.9455 0.1727
Outcomes 0.09096 0.04548     2 40.083  0.2647 0.7688

# Extract p-values
p_time <- a_m2["time", "Pr(>F)"]
p_outcomes <- a_m2["Outcomes", "Pr(>F)"]

annot_label <- paste0(
  "p[time] == ", format(round(p_time, 3), nsmall = 3),
  "*';'*~p[outcomes] == ", format(round(p_outcomes, 3), nsmall = 3)
)

# ==========================
# Diagnostics
# ==========================
plot(mod1_m2)

plot(fitted(mod1_m2), residuals(mod1_m2),
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs Fitted")

qqnorm(residuals(mod1_m2))
qqline(residuals(mod1_m2), col = "blue")

acf(residuals(mod1_m2), main = "ACF of residuals")

qqnorm(ranef(mod1_m2)$ID[, "(Intercept)"], main = "QQ plot of random intercepts")
qqline(ranef(mod1_m2)$ID[, "(Intercept)"], col = "blue")

# ==========================
# Final Plot
# ==========================
Figure_1m_2 <- ggplot(F1m_2 %>% filter(value > 0), aes(
  x = ADR2_time, y = value, fill = Outcomes,
  group = interaction(ADR2_time, Outcomes)
)) +
  geom_boxplot(position = position_dodge(0.7), width = 0.6) +
  scale_y_log10() +
  labs(
    title = "Association between ADRB2 expression in monocytes and outcomes",
    x = "Time of measurement",
    y = "Percentage of monocytes expressing ADRB2\n(log10 scale)"
  ) +
  geom_text(data = count_data_m2, aes(
    x = ADR2_time,
    y = min(F1m_2$value[F1m_2$value > 0], na.rm = TRUE) * 0.5,
    label = paste0("n=", n),
    group = interaction(ADR2_time, Outcomes)
  ),
  position = position_dodge(0.7),
  size = 4,
  vjust = 1) +
  scale_x_discrete(labels = c(
    "ADR2m_J0" = "implantation",
    "ADR2m_J3_J5" = "day 3 to 5",
    "ADR2m_JS" = "explantation"
  )) +
    scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",                   
    "Bridge to transplant or LVAD" = "#377EB8",   
    "Death" = "#F8766D"                           
  )) +
  annotate(
    "text",
    x = 1,
    y = max(F1m_2$value, na.rm = TRUE) * 1.2,
    label = annot_label,
    parse = TRUE,
    size = 4,
    hjust = 0
  )

Fig_1m <- Figure_1m # + theme(plot.title = element_blank())
Fig_1m_2 <- Figure_1m_2 # + theme(plot.title = element_blank())

Figure1m_tot <- Fig_1m + Fig_1m_2 +
  plot_annotation(title = "Association between percentage of monocytes expressing ADRB and Outcome")

Figure1m_tot

Statistical Analysis

To explore the distribution of ADRB1 and ADRB2 expression over time in relation to clinical outcomes, boxplots were generated at three timepoints (implantation, day 3–5, and explantation). To evaluate the temporal dynamics of ADRB1 and ADRB2 expression across clinical outcomes, we used linear mixed-effects models with time, outcome group (ECMO weaning, bridge to LVAD/transplant, or death), and their interaction as fixed effects. Random intercepts and slopes for time were included to account for intra-patient variability. When appropriate, receptor expression values were log10-transformed. Interaction terms were tested using likelihood ratio tests between full and reduced models, and Type III ANOVA with Kenward-Roger approximation was used for fixed effects.

Results (T4 lymphocytes)

In patients who were successfully weaned from ECMO, both ADRB1 and ADRB2 expression in T4 lymphocytes showed an increase from implantation to explantation. In contrast, patients who died during ECMO exhibited lower or decreasing receptor expression over time. Although these trends were visually distinguishable, the interaction between time and outcome group was not statistically significant (ADRB1: p = 0.58; ADRB2: p = 0.76). Nevertheless, we observed a significant main effect of time for ADRB1 (p = 0.02) and a trend for ADRB2 (p = 0.07), alongside borderline differences across outcome groups (ADRB1: p = 0.07; ADRB2: p = 0.05).

Results (monocytes)

In monocytes, ADRB1 and ADRB2 expression did not show any consistent time-dependent pattern across outcome groups. Boxplots revealed no evident increase or decrease over time, regardless of clinical trajecotry. These observations were supported by linear mixed-effects models, which showed no significant interaction between time and outcome group (ADRB1: pinteraction = 0.64; ADRB2: pinteraction = 0.76). Furthermore, neither time nor outcome group had significant main effects (ADRB1: p[time] = 0.95, p[outcomes] = 0.65; ADRB2: p[time] = 0.94, p[outcomes] = 0.57).

Figure 2: Association between VTI and ADRB in Lymphocyte at day of weaning

#ADRB1
F2_1_a <- df %>%
select(`jsev_epr_sev_t1_1,5l_itv`, ADR1_JS)%>%
  mutate(
    VTI_JS = case_when(
      `jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
      `jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
    )
  )%>%
  na.omit()
   

 # Mann-Whitney test for VTI and ADRB1
mann_whitney_VTI <- wilcox.test(ADR1_JS ~ VTI_JS, data = F2_1_a)

# Extract and format the p-value
p_value_VTI <- sprintf("%.3f", mann_whitney_VTI$p.value)

count_data <- F2_1_a %>%
  group_by(VTI_JS) %>%
  summarise(n = n())

# Boxplot for ADRB1 with p-value annotation 
Figure2a <- ggplot(F2_1_a, aes(x = VTI_JS, y = log10(ADR1_JS), fill = VTI_JS)) +
  geom_boxplot() +
  geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
  geom_text(data = count_data,
            aes(x = VTI_JS, y = min(log10(F2_1_a$ADR1_JS), na.rm = TRUE) * 1.25,
                label = paste0("n=", n)),
            size = 4,
            inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
  labs(
    title = "A",
    x = "VTI measured at 1.5 L/min",
    y = "ADRB1 in lymphocyte measured at ECMO weaning \n (log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text",
           x = 1.5,
           y = max(log10(F2_1_a$ADR1_JS), na.rm = TRUE),
           label = paste("p:", p_value_VTI),
           size = 5, color = "black", hjust = 0.5)


#ADRB2
F2_2_a <- df %>%
select(`jsev_epr_sev_t1_1,5l_itv`, ADR2_JS)%>%
  mutate(
    VTI_JS = case_when(
      `jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
      `jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
    )
  )%>%
  na.omit()
   

 # Mann-Whitney test for VTI and ADRB2
mann_whitney_VTI_2 <- wilcox.test(log10(ADR2_JS) ~ VTI_JS, data = F2_2_a)

# Extract and format the p-value
p_value_VTI_2 <- sprintf("%.3f", mann_whitney_VTI_2$p.value)

count_data <- F2_2_a %>%
  group_by(VTI_JS) %>%
  summarise(n = n(), .groups = "drop")

# Boxplot for ADRB2 with p-value annotation 
Figure2b <- ggplot(F2_2_a, aes(x = VTI_JS, y = log10(ADR2_JS), fill = VTI_JS)) +
  geom_boxplot() +
  geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
  scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
  labs(
    title = "B",
    x = "VTI measured at 1.5L /min",
    y = "ADRB2 in lymphocyte measured at ECMO weaning \n (log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate(
    "text",
    x = 1.5,
    y = max(log10(F2_2_a$ADR2_JS), na.rm = TRUE),
    label = paste("p:", p_value_VTI_2),
    size = 5,
    color = "black",
    hjust = 0.5
  ) +
  geom_text(
    data = count_data,
    aes(
      x = VTI_JS,
      y = min(log10(F2_2_a$ADR2_JS), na.rm = TRUE) * 1.25,  
      label = paste0("n=", n)
    ),
    inherit.aes = FALSE,
    size = 4
  )

Figure2_tot <- Figure2a | Figure2b +
  plot_annotation(
    title = "Association between VTI at weaning and ADRB1/ADRB2 expression in lymphocyte",
    tag_levels = c("A","B")  
  )

Figure2_tot

Figure 2_m: Association between VTI and ADRB in Monocytes at day of weaning

# ADRB1

F2_1_m <- df %>%
  select(`jsev_epr_sev_t1_1,5l_itv`, ADR1m_JS) %>%
  mutate(
    VTI_JS = case_when(
      `jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
      `jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
    ),
    VTI_JS = factor(VTI_JS, levels = c("≤12 cm", ">12 cm"))
  ) %>%
  na.omit()

# Mann–Whitney test
mw_test_adr1m <- wilcox.test(log10(ADR1m_JS) ~ VTI_JS, data = F2_1_m)

# Extract p-value
p_value_adr1m <- sprintf("%.3f", mw_test_adr1m$p.value)

count_data_m <- F2_1_m %>%
  group_by(VTI_JS) %>%
  summarise(n = n())

# Plot
Figure2m_A <- ggplot(F2_1_m, aes(x = VTI_JS, y = log10(ADR1m_JS), fill = VTI_JS)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_m,
            aes(x = VTI_JS,
                y = -2.5,
                label = paste0("n=", n)),
            size = 4,
            inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
  labs(
    title = "A",
    x = "VTI measured at 1.5L /min",
    y = "ADRB1 in monocyte at ECMO weaning\n(log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate(
    "text",
    x = 1.5,
    y = max(log10(F2_1_m$ADR1m_JS), na.rm = TRUE) * 1.05,
    label = paste("p:", p_value_adr1m),
    size = 5
  ) +
  coord_cartesian(clip = "off")

# ADRB2

F2_2_m <- df %>%
  select(`jsev_epr_sev_t1_1,5l_itv`, ADR2m_JS) %>%
  mutate(
    VTI_JS = case_when(
      `jsev_epr_sev_t1_1,5l_itv` > 12 ~ ">12 cm",
      `jsev_epr_sev_t1_1,5l_itv` <= 12 ~ "≤12 cm"
    ),
    VTI_JS = factor(VTI_JS, levels = c("≤12 cm", ">12 cm"))
  ) %>%
  na.omit()

# Mann–Whitney test
mw_test_adr2m <- wilcox.test(log10(ADR2m_JS) ~ VTI_JS, data = F2_2_m)

# Extract p-value
p_value_adr2m <- sprintf("%.3f", mw_test_adr2m$p.value)

count_data_m2 <- F2_2_m %>%
  group_by(VTI_JS) %>%
  summarise(n = n())

# Plot
Figure2m_B <- ggplot(F2_2_m, aes(x = VTI_JS, y = log10(ADR2m_JS), fill = VTI_JS)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_m2,
            aes(x = VTI_JS,
                y = min(log10(F2_2_m$ADR2m_JS), na.rm = TRUE) * 1.25,
                label = paste0("n=", n)),
            size = 4,
            inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤12 cm" = "#0072B2", ">12 cm" = "#D55E00")) +
  labs(
    title = "B",
    x = "VTI measured at 1.5L /min",
    y = "ADRB2 in monocyte at ECMO weaning\n(log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate(
    "text",
    x = 1.5,
    y = max(log10(F2_2_m$ADR2m_JS), na.rm = TRUE) * 1.05,
    label = paste("p:", p_value_adr2m),
    size = 5
  )

# Combined plot

Figure2m_tot <- Figure2m_A | Figure2m_B +
  plot_annotation(
    title = "Association between VTI at weaning and ADRB1/ADRB2 expression in monocyte",
    tag_levels = "A"
  )

Figure2m_tot

Statistical Analysis

To assess the association between ADRB1 and ADRB2 expression at the time of ECMO weaning and left ventricular outflow tract velocity-time integral (VTI) measured at 1.5 L/min ECMO flow, patients were stratified according to a VTI threshold of 12 cm. Mann–Whitney U tests were used to compare log-transformed ADRB expression between the two groups (≤12 cm vs >12 cm). Results were visualized using boxplots.

Results

In lymphocytes, both ADRB1 and ADRB2 expression levels at the time of ECMO weaning were significantly lower in patients with a suboptimal VTI (≤12 cm), a threshold commonly used to assess readiness for explantation. This difference was statistically significant for both markers (ADRB1: p = 0.04; ADRB2: p = 0.04), suggesting a potential association between impaired hemodynamic recovery and reduced β-adrenergic receptor expression.

In contrast, no statistically significant differences in ADRB1 or ADRB2 expression in monocytes were observed between patients with VTI ≤12 cm and those with VTI >12 cm (ADRB1: p = 0.73; ADRB2: p = 0.79).

ADR and inotropics

ADR ~ dobutamin

### Correlation between ADRB1 in Lymphocyte and total dobutamine dose

F3c_a <- df %>%
  select(ADR1_JS, dose_tot_dobu) %>%
  filter(!is.na(ADR1_JS), !is.na(dose_tot_dobu))

cor_test_adr1 <- cor.test(F3c_a$ADR1_JS, F3c_a$dose_tot_dobu, method = "spearman")
rho_adr1 <- sprintf("%.4f", cor_test_adr1$estimate)
pval_adr1 <- sprintf("p = %.3f", cor_test_adr1$p.value)

Fig_3c <- ggplot(F3c_a, aes(x = dose_tot_dobu, y = ADR1_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(title = "Correlation between dobutamine dose and ADRB1 (lymphocytes)",
       x = "Dobutamine dose (mg)",
       y = "ADRB1 expression on T helper lymphocytes") +
  annotate("text", 
           x = min(F3c_a$dose_tot_dobu, na.rm = TRUE), 
           y = max(F3c_a$ADR1_JS, na.rm = TRUE), 
           label = paste("Spearman rho:", rho_adr1, "\n", pval_adr1),
           hjust = 0, size = 5, color = "black") +
  theme_minimal()


### Correlation between ADRB2 in Lymphocyte and total dobutamine dose

F3c_b <- df %>%
  select(ADR2_JS, dose_tot_dobu) %>%
  filter(!is.na(ADR2_JS), !is.na(dose_tot_dobu))

cor_test_adr2 <- cor.test(F3c_b$ADR2_JS, F3c_b$dose_tot_dobu, method = "spearman")
rho_adr2 <- sprintf("%.4f", cor_test_adr2$estimate)
pval_adr2 <- sprintf("p = %.3f", cor_test_adr2$p.value)

Fig_3d <- ggplot(F3c_b, aes(x = dose_tot_dobu, y = ADR2_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(title = "Correlation between dobutamine dose and ADRB2 (lymphocytes)",
       x = "Dobutamine dose (mg)",
       y = "ADRB2 expression on T helper lymphocytes") +
  annotate("text", 
           x = min(F3c_b$dose_tot_dobu, na.rm = TRUE), 
           y = max(F3c_b$ADR2_JS, na.rm = TRUE), 
           label = paste("Spearman rho:", rho_adr2, "\n", pval_adr2),
           hjust = 0, size = 5, color = "black") +
  theme_minimal()

### Correlation between ADRB1 in Monocyte and total dobutamine dose

F3c_m1 <- df %>%
  select(ADR1m_JS, dose_tot_dobu) %>%
  filter(!is.na(ADR1m_JS), !is.na(dose_tot_dobu))

cor_test_adr1m <- cor.test(F3c_m1$ADR1m_JS, F3c_m1$dose_tot_dobu, method = "spearman")
rho_adr1m <- sprintf("%.4f", cor_test_adr1m$estimate)
pval_adr1m <- sprintf("p = %.3f", cor_test_adr1m$p.value)

Fig_3c_m <- ggplot(F3c_m1, aes(x = dose_tot_dobu, y = ADR1m_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(title = "Correlation between dobutamine dose and ADRB1 (monocytes)",
       x = "Dobutamine dose (mg)",
       y = "ADRB1 expression on monocytes at ECMO weaning") +
  annotate("text", 
           x = min(F3c_m1$dose_tot_dobu, na.rm = TRUE), 
           y = max(F3c_m1$ADR1m_JS, na.rm = TRUE), 
           label = paste("Spearman rho:", rho_adr1m, "\n", pval_adr1m),
           hjust = 0, size = 5, color = "black") +
  theme_minimal()

### Correlation between ADRB2 in Monocyte and total dobutamine dose

F3c_m2 <- df %>%
  select(ADR2m_JS, dose_tot_dobu) %>%
  filter(!is.na(ADR2m_JS), !is.na(dose_tot_dobu))

cor_test_adr2m <- cor.test(F3c_m2$ADR2m_JS, F3c_m2$dose_tot_dobu, method = "spearman")
rho_adr2m <- sprintf("%.4f", cor_test_adr2m$estimate)
pval_adr2m <- sprintf("p = %.3f", cor_test_adr2m$p.value)

Fig_3d_m <- ggplot(F3c_m2, aes(x = dose_tot_dobu, y = ADR2m_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(title = "Correlation between dobutamine dose and ADRB2 (monocytes)",
       x = "Dobutamine dose (mg)",
       y = "ADRB2 expression on monocytes at ECMO weaning") +
  annotate("text", 
           x = min(F3c_m2$dose_tot_dobu, na.rm = TRUE), 
           y = max(F3c_m2$ADR2m_JS, na.rm = TRUE), 
           label = paste("Spearman rho:", rho_adr2m, "\n", pval_adr2m),
           hjust = 0, size = 5, color = "black") +
  theme_minimal()

Statistical analysis

The association between cumulative dobutamine dose and the percentage of T-helper lymphocytes/monocyte expressing ADRB was assessed using Spearman’s rank correlation coefficient.

Results

A weak negative correlation was observed between ADRB expression in lymphocytes and total dobutamine dose, which did not reach statistical significance. In contrast, both ADRB1 and ADRB2 expression in monocytes at the time of ECMO weaning were moderately correlated with cumulative dobutamine dose (Spearman ρ = 0.52, p = 0.021 for ADRB1; ρ = 0.54, p = 0.015 for ADRB2), indicating a significant positive association.

Figure 3: Association of cumulative dobutamine and ADRB in Lymphocyte at day of weaning

# ADRB1

F3_a <- df %>%
  select(dose_tot_dobu, ADR1_JS) %>%
  mutate(dobu_cut = case_when(
    dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
    dose_tot_dobu <= median(dose_tot_dobu, na.rm = TRUE) ~ "≤ median"
  )) %>%
  filter(!is.na(ADR1_JS), !is.na(dose_tot_dobu))%>%
  mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))

# Test di Mann–Whitney
mw_adr1 <- wilcox.test(log10(ADR1_JS) ~ dobu_cut, data = F3_a)
p_val_adr1 <- sprintf("%.3f", mw_adr1$p.value)

count_data_a <- F3_a %>%
  group_by(dobu_cut) %>%
  summarise(n = n(), .groups = "drop")

# Figure ADR1
Figure3a <- ggplot(F3_a, aes(x = dobu_cut, y = log10(ADR1_JS), fill = dobu_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
  geom_text(data = count_data_a,
            aes(x = dobu_cut, y = min(log10(F3_a$ADR1_JS), na.rm = TRUE) * 1.1,
                label = paste0("n=", n)),
            size = 4,
            inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
  labs(
    title = "A",
    x = "Cumulative Dobutamine (cut at median)",
    y = "ADRB1 at ECMO weaning (log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text",
           x = 1.5,
           y = max(log10(F3_a$ADR1_JS), na.rm = TRUE) * 1.15,
           label = paste("p:", p_val_adr1),
           size = 5)

# ADRB2

F3_b <- df %>%
  select(dose_tot_dobu, ADR2_JS) %>%
  mutate(dobu_cut = case_when(
    dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
    dose_tot_dobu <= median(dose_tot_dobu, na.rm = TRUE) ~ "≤ median"
  )) %>%
  filter(!is.na(ADR2_JS), !is.na(dose_tot_dobu))%>%
  mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))

# Test
mw_adr2 <- wilcox.test(log10(ADR2_JS) ~ dobu_cut, data = F3_b)
p_val_adr2 <- sprintf("%.3f", mw_adr2$p.value)

count_data_b <- F3_b %>%
  group_by(dobu_cut) %>%
  summarise(n = n(), .groups = "drop")

# Figure ADR2
Figure3b <- ggplot(F3_b, aes(x = dobu_cut, y = log10(ADR2_JS), fill = dobu_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
  geom_text(data = count_data_b,
            aes(x = dobu_cut, y = min(log10(F3_b$ADR2_JS), na.rm = TRUE) * 1.1,
                label = paste0("n=", n)),
            size = 4,
            inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
  labs(
    title = "B",
    x = "Cumulative Dobutamine (cut at median)",
    y = "ADRB2 at ECMO weaning (log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text",
           x = 1.5,
           y = max(log10(F3_b$ADR2_JS), na.rm = TRUE) * 1.15,
           label = paste("p:", p_val_adr2),
           size = 5)

# === Combine ===
Figure3_dobu_cum <- Figure3a + Figure3b +
  plot_annotation(title = "Association between cumulative dobutamine and ADR expression in lymphocytes")

Figure3_dobu_cum

Figure 3_m: Association of cumulative dobutamine and ADRB in Monocyte at day of weaning

## ADRB1

# Prepare data
F3_m_a <- df %>%
  select(dose_tot_dobu, ADR1m_JS) %>%
  mutate(
    dobu_cut = case_when(
      dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
      TRUE ~ "≤ median"
    )
  ) %>%
  filter(!is.na(ADR1m_JS), !is.na(dose_tot_dobu)) %>%
  mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))


# Mann–Whitney test
mw_test_adr1m <- wilcox.test(log10(ADR1m_JS) ~ dobu_cut, data = F3_m_a)
p_value_adr1m <- sprintf("%.3f", mw_test_adr1m$p.value)

# Count per group
count_data_adr1m <- F3_m_a %>%
  group_by(dobu_cut) %>%
  summarise(n = n())

# Plot
Figure3_dobu_cum_a <- ggplot(F3_m_a, aes(x = dobu_cut, y = log10(ADR1m_JS), fill = dobu_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
  geom_text(
    data = F3_m_a %>% group_by(dobu_cut) %>% summarise(n = n()) %>%
      mutate(y_n_label = -2.5),  # y personalizzata direttamente qui
    aes(x = dobu_cut, y = y_n_label, label = paste0("n=", n)),
    inherit.aes = FALSE, size = 4
  ) +
  scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
  labs(title = "A", x = "Cumulative Dobutamine", y = "ADRB1 in monocytes at weaning (log10)") +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text",
           x = 1.5,
           y = max(log10(F3_m_a$ADR1m_JS), na.rm = TRUE),
           label = paste("p:", p_value_adr1m),
           size = 5)


##ADRB2

# Prepare data
F3_m_b <- df %>%
  select(dose_tot_dobu, ADR2m_JS) %>%
  mutate(
    dobu_cut = case_when(
      dose_tot_dobu > median(dose_tot_dobu, na.rm = TRUE) ~ "> median",
      TRUE ~ "≤ median"
    )
  ) %>%
  filter(!is.na(ADR2m_JS), !is.na(dose_tot_dobu))%>%
  mutate(dobu_cut = factor(dobu_cut, levels = c("≤ median", "> median")))

# Mann–Whitney test
mw_test_adr2m <- wilcox.test(log10(ADR2m_JS) ~ dobu_cut, data = F3_m_b)
p_value_adr2m <- sprintf("%.3f", mw_test_adr2m$p.value)

# Count per group
count_data_adr2m <- F3_m_b %>%
  group_by(dobu_cut) %>%
  summarise(n = n())

# Plot
Figure3_dobu_cum_b <- ggplot(F3_m_b, aes(x = dobu_cut, y = log10(ADR2m_JS), fill = dobu_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.7, size = 2) +
  geom_text(data = count_data_adr2m, aes(x = dobu_cut, y = min(log10(F3_m_b$ADR2m_JS), na.rm = TRUE) * 1.25,
                                         label = paste0("n=", n)), inherit.aes = FALSE, size = 4) +
  scale_fill_manual(values = c("≤ median" = "#009E73", "> median" = "#CC79A7")) +
  labs(title = "B", x = "Cumulative Dobutamine", y = "ADRB2 in monocytes at weaning (log10)") +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text", x = 1.5, y = max(log10(F3_m_b$ADR2m_JS), na.rm = TRUE), label = paste("p:", p_value_adr2m), size = 5)

# Cumulative
Figure3m_dobu_cum <- Figure3_dobu_cum_a + Figure3_dobu_cum_b +
  plot_annotation(title = "Association between cumulative dobutamine and ADR expression in monocytes")

Figure3m_dobu_cum

Statistical Analysis

As evolution of some patients allowed an early weaning before D5, we considered that D3-D5 and D weaning are the same time point. To explore whether cumulative exposure to inotropic stimulation was associated with β-adrenergic receptor modulation over time, we investigated the relationship between the total cumulative dose of dobutamine (expressed in mg) and ADRB1/ADRB2 expression at ECMO weaning day. This timepoint was selected as it reflects the net effect of continuous pharmacological stimulation throughout the ECMO course, thereby allowing the assessment of potential receptor adaptation (e.g., up- or downregulation).

First, we performed a correlation analysis between cumulative dobutamine dose and ADRB1/ADRB2 expression using Spearman’s rank correlation coefficient, due to the non-normal distribution of the data. Second, patients were stratified according to the median of cumulative dobutamine dose. Expression levels of ADRB1 and ADRB2 (in both CD4+ lymphocytes and monocytes) at ECMO weaning were compared between patients above versus below the median using the Mann–Whitney U test. A two-sided p-value < 0.05 was considered statistically significant.

Results

Among lymphocytes, there was no significant difference in ADRB expression between patients with higher versus lower cumulative dobutamine exposure. ADRB1 levels were comparable across groups, while a trend toward lower ADRB2 expression was observed in patients who received higher doses of dobutamine.

In contrast, monocyte expression patterns followed a more consistent trend. Both ADRB1 and ADRB2 expression levels tended to be higher in patients exposed to higher cumulative doses of dobutamine, although this association did not reach statistical significance.

ADR ~ noradrenalin

### Correlation between ADRB1 in Lymphocyte and total noradrenalin dose

df4 <- df %>%
  dplyr::select(dose_tot_NAd, ADR1_JS, Weight) %>%
  filter(!is.na(ADR1_JS), !is.na(dose_tot_NAd)) %>%
  mutate(NAD_per_kg = dose_tot_NAd / Weight)

# Correlation
correlation_adr1 <- sprintf("%.4f", cor(df4$ADR1_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))

# Plot ADRB1
Figure_4_NADcum_A <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR1_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Correlation between noradrenaline dose and ADRB1 (lymphocytes)",
    x = "Cumulative noradrenaline dose (mg)",
    y = "ADRB1 expression at ECMO weaning (JS)"
  ) +
  annotate(
    "text",
    x = min(df4$dose_tot_NAd, na.rm = TRUE),
    y = max(df4$ADR1_JS, na.rm = TRUE) * 0.95,
    label = paste("Spearman correlation:", correlation_adr1),
    hjust = 0,
    size = 5,
    color = "black"
  ) +
  theme_minimal()


print(Figure_4_NADcum_A)

### Correlation between ADRB2 in Lymphocyte and total noradrenalin dose

df4 <- df %>%
  dplyr::select(dose_tot_NAd, ADR2_JS, Weight) %>%
  filter(!is.na(ADR2_JS), !is.na(dose_tot_NAd)) %>%
  mutate(NAD_per_kg = dose_tot_NAd / Weight)

correlation_adr2 <- sprintf("%.4f", cor(df4$ADR2_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))

Figure_4_NADcum_B <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR2_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Correlation between noradrenaline dose and ADRB2 (lymphocytes)",
    x = "Cumulative noradrenaline dose (mg)",
    y = "ADRB2 expression at ECMO weaning (JS)"
  ) +
  annotate(
    "text",
    x = min(df4$dose_tot_NAd, na.rm = TRUE),
    y = max(df4$ADR2_JS, na.rm = TRUE) * 0.95,
    label = paste("Spearman correlation:", correlation_adr2),
    hjust = 0,
    size = 5,
    color = "black"
  ) +
  theme_minimal()

print(Figure_4_NADcum_B)

### Correlation between ADRB1 in Monocyte and total noradrenalin dose

df4 <- df %>%
  select(ADR1m_JS, dose_tot_NAd, Weight) %>%
  filter(!is.na(ADR1m_JS), !is.na(dose_tot_NAd)) %>%
  mutate(NAD_per_kg = dose_tot_NAd / Weight)

# Spearman correlation
cor_adr1m <- sprintf("%.4f", cor(df4$ADR1m_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))

# Plot
Figure4_NADmono_A <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR1m_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Correlation between Noradrenaline dose and ADRB1 (monocytes)",
    x = "Cumulative noradrenaline dose (mg)",
    y = "Monocyte ADRB1 expression at ECMO weaning"
  ) +
  annotate("text",
           x = min(df4$dose_tot_NAd, na.rm = TRUE),
           y = max(df4$ADR1m_JS, na.rm = TRUE) * 0.95,
           label = paste("Spearman correlation:", cor_adr1m),
           hjust = 0, size = 5, color = "black") +
  theme_minimal()

print(Figure4_NADmono_A)

### Correlation between ADRB2 in Monocyte and total noradrenalin dose

df4 <- df %>%
  select(ADR2m_JS, dose_tot_NAd, Weight) %>%
  filter(!is.na(ADR2m_JS), !is.na(dose_tot_NAd)) %>%
  mutate(NAD_per_kg = dose_tot_NAd / Weight)

# Spearman correlation
cor_adr2m <- sprintf("%.4f", cor(df4$ADR2m_JS, df4$dose_tot_NAd, method = "spearman", use = "complete.obs"))

# Plot
Figure4_NADmono_B <- ggplot(df4, aes(x = dose_tot_NAd, y = ADR2m_JS)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = TRUE) +
  labs(
    title = "Correlation between Noradrenaline dose and ADRB2 (monocytes)",
    x = "Cumulative noradrenaline dose (mg)",
    y = "Monocyte ADRB2 expression at ECMO weaning"
  ) +
  annotate("text",
           x = min(df4$dose_tot_NAd, na.rm = TRUE),
           y = max(df4$ADR2m_JS, na.rm = TRUE) * 0.95,
           label = paste("Spearman correlation:", cor_adr2m),
           hjust = 0, size = 5, color = "black") +
  theme_minimal()

print(Figure4_NADmono_B)

Results

No significant correlations were observed between cumulative noradrenaline dose and ADRB1 or ADRB2 expression in either T helper lymphocytes (r = –0.09 and –0.29, respectively) or monocytes (r = –0.18 and 0.02, respectively).

Figure 4: Association of cumulative noradrenaline and ADRB in Lymphocyte at day of weaning

### ADRB1 in T lymphocytes ~ cumulative NAd (median)

F5a <- df %>%
  mutate(
    ADR1_JS_clean = case_when(
      !is.na(ADR1_JS) ~ ADR1_JS,
      is.na(ADR1_JS) & !is.na(ADR1_J3_J5) ~ ADR1_J3_J5,
      TRUE ~ NA_real_
    ),
    NAd_group = case_when(
      dose_tot_NAd > median(dose_tot_NAd, na.rm = TRUE) ~ "> median",
      dose_tot_NAd <= median(dose_tot_NAd, na.rm = TRUE) ~ "≤ median",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(ADR1_JS_clean), !is.na(NAd_group)) %>%
  mutate(NAd_group = factor(NAd_group, levels = c("≤ median", "> median")))

# p-value (Mann–Whitney)
wilcox_test_adr1 <- wilcox.test(ADR1_JS_clean ~ NAd_group, data = F5a)
p_value_adr1 <- sprintf("%.3f", wilcox_test_adr1$p.value)

count_data_adr1 <- F5a %>%
  group_by(NAd_group) %>%
  summarise(n = n())

# Boxplot
Figure_5a <- ggplot(F5a, aes(x = NAd_group, y = ADR1_JS_clean, fill = NAd_group)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(
    data = count_data_adr1,
    aes(x = NAd_group, y = min(F5a$ADR1_JS_clean, na.rm = TRUE) * 1.2,
        label = paste0("n=", n)),
    size = 4,
    inherit.aes = FALSE
  ) +
  scale_fill_manual(values = c("≤ median" = "#66C2AC", "> median" = "#D73027")) +
  labs(
    title = "A",
    x = "Cumulative noradrenaline dose",
    y = "ADRB1 expression at ECMO weaning"
  ) +
  annotate("text",
           x = 1.5,
           y = max(F5a$ADR1_JS_clean, na.rm = TRUE) * 1.05,
           label = paste("p:", p_value_adr1),
           size = 5, color = "black") +
  theme_minimal() +
  theme(legend.position = "none")


### ADRB2 in T lymphocytes ~ cumulative NAd (median)

F5b <- df %>%
  mutate(
    ADR2_JS_clean = case_when(
      !is.na(ADR2_JS) ~ ADR2_JS,
      is.na(ADR2_JS) & !is.na(ADR2_J3_J5) ~ ADR2_J3_J5,
      TRUE ~ NA_real_
    ),
    NAd_group = case_when(
      dose_tot_NAd > median(dose_tot_NAd, na.rm = TRUE) ~ "> median",
      dose_tot_NAd <= median(dose_tot_NAd, na.rm = TRUE) ~ "≤ median",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(ADR2_JS_clean), !is.na(NAd_group)) %>%
  mutate(NAd_group = factor(NAd_group, levels = c("≤ median", "> median")))

# Test di Mann–Whitney
wilcox_test_adr2 <- wilcox.test(ADR2_JS_clean ~ NAd_group, data = F5b)
p_value_adr2 <- sprintf("%.3f", wilcox_test_adr2$p.value)

# Conta per gruppo
count_data_adr2 <- F5b %>%
  group_by(NAd_group) %>%
  summarise(n = n())

# Boxplot
Figure_5b <- ggplot(F5b, aes(x = NAd_group, y = ADR2_JS_clean, fill = NAd_group)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(
    data = count_data_adr2,
    aes(x = NAd_group, y = min(F5b$ADR2_JS_clean, na.rm = TRUE) * 1.2,
        label = paste0("n=", n)),
    size = 4,
    inherit.aes = FALSE
  ) +
  scale_fill_manual(values = c("≤ median" = "#66C2AC", "> median" = "#D73027")) +
  labs(
    title = "B",
    x = "Cumulative noradrenaline dose",
    y = "ADRB2 expression at ECMO weaning"
  ) +
  annotate("text",
           x = 1.5,
           y = max(F5b$ADR2_JS_clean, na.rm = TRUE) * 1.05,
           label = paste("p:", p_value_adr2),
           size = 5, color = "black") +
  theme_minimal() +
  theme(legend.position = "none")


Figure_5a <- Figure_5a + theme(plot.title = element_blank())
Figure_5b <- Figure_5b + theme(plot.title = element_blank())

# Combinazione dei due plot
Figure4_tot <- Figure_5a + Figure_5b +
  plot_annotation(
    title = "Association between cumulative noradrenaline dose and β-adrenergic receptor expression on T helper lymphocytes"
  )

Figure4_tot

Figure 4m: Association of cumulative noradrenaline and ADRB in Monocyte at day of weaning

cut_median_nad <- median(df$dose_tot_NAd, na.rm = TRUE)

#ADRB1m
F4_m_a <- df %>%
  select(dose_tot_NAd, ADR1m_JS, ADR1m_J3_J5) %>%
  mutate(
    ADR1m_final = case_when(
      !is.na(ADR1m_JS) ~ ADR1m_JS,
      is.na(ADR1m_JS) & !is.na(ADR1m_J3_J5) ~ ADR1m_J3_J5,
      TRUE ~ NA_real_
    ),
    NAd_cut = case_when(
      dose_tot_NAd <= cut_median_nad ~ "≤ median",
      dose_tot_NAd > cut_median_nad ~ "> median",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(ADR1m_final), !is.na(NAd_cut)) %>%
  mutate(NAd_cut = factor(NAd_cut, levels = c("≤ median", "> median")))

# P-value
p_val_adr1m <- wilcox.test(ADR1m_final ~ NAd_cut, data = F4_m_a)$p.value
p_text_adr1m <- paste("p:", format(round(p_val_adr1m, 3), nsmall = 3))

# Count
count_data_m1 <- F4_m_a %>%
  group_by(NAd_cut) %>%
  summarise(n = n())

# Plot ADRB1m
Figure4m_A <- ggplot(F4_m_a, aes(x = NAd_cut, y = log10(ADR1m_final), fill = NAd_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_m1,
            aes(x = NAd_cut, y = min(log10(F4_m_a$ADR1m_final), na.rm = TRUE) * 1.2,
                label = paste0("n=", n)),
            size = 4, inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤ median" = "#66C2A5", "> median" = "#D73027")) +
  labs(
    title = "A",
    x = "Cumulative norepinephrine dose",
    y = "Monocyte ADRB1 expression at ECMO weaning\n(log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text", x = 1.5,
           y = max(log10(F4_m_a$ADR1m_final), na.rm = TRUE) * 1.15,
           label = p_text_adr1m,
           size = 5)


# ADRB2m

F4_m_b <- df %>%
  select(dose_tot_NAd, ADR2m_JS, ADR2m_J3_J5) %>%
  mutate(
    ADR2m_final = case_when(
      !is.na(ADR2m_JS) ~ ADR2m_JS,
      is.na(ADR2m_JS) & !is.na(ADR2m_J3_J5) ~ ADR2m_J3_J5,
      TRUE ~ NA_real_
    ),
    NAd_cut = case_when(
      dose_tot_NAd <= cut_median_nad ~ "≤ median",
      dose_tot_NAd > cut_median_nad ~ "> median",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(ADR2m_final), !is.na(NAd_cut)) %>%
  mutate(NAd_cut = factor(NAd_cut, levels = c("≤ median", "> median")))

# P-value
p_val_adr2m <- wilcox.test(ADR2m_final ~ NAd_cut, data = F4_m_b)$p.value
p_text_adr2m <- paste("p:", format(round(p_val_adr2m, 3), nsmall = 3))

# Count
count_data_m2 <- F4_m_b %>%
  group_by(NAd_cut) %>%
  summarise(n = n())

# Plot ADRB2m
Figure4m_B <- ggplot(F4_m_b, aes(x = NAd_cut, y = log10(ADR2m_final), fill = NAd_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_m2,
            aes(x = NAd_cut, y = min(log10(F4_m_b$ADR2m_final), na.rm = TRUE) * 1.2,
                label = paste0("n=", n)),
            size = 4, inherit.aes = FALSE) +
  scale_fill_manual(values = c("≤ median" = "#66C2AC", "> median" = "#D73027")) +
  labs(
    title = "B",
    x = "Cumulative norepinephrine dose",
    y = "Monocyte ADRB2 expression at ECMO weaning\n(log10 scale)"
  ) +
  theme_minimal() +
  theme(legend.position = "none") +
  annotate("text", x = 1.5,
           y = max(log10(F4_m_b$ADR2m_final), na.rm = TRUE) * 1.15,
           label = p_text_adr2m,
           size = 5)

# Combine plots

Figure4_mono <- Figure4m_A + Figure4m_B +
  plot_annotation(title = "Association between cumulative norepinephrine dose and ADRB1/2 expression in monocytes at ECMO weaning")

Figure4_mono

Results

Among lymphocytes, both ADRB1 and ADRB2 expression at ECMO weaning tended to be lower in patients receiving a cumulative noradrenaline dose above the median. However, these differences did not reach statistical significance.

Among monocytes, ADRB1 expression at ECMO weaning was significantly lower in patients receiving a cumulative norepinephrine dose above the median compared to those below (p = 0.039). A similar trend was observed for ADRB2 expression, although the difference did not reach statistical significance (p = 0.077).

While Spearman correlation analyses did not reveal a significant association between cumulative norepinephrine dose and ADRB1/ADRB2 expression in monocytes, patients with doses above the median exhibited significantly lower ADRB1 levels at ECMO weaning (p = 0.039). This suggests that potential receptor downregulation may occur only above a certain threshold of catecholaminergic stimulation.

ADRB ~ VIS Score at day 0 in Lymphocytes

df <- df %>%
  mutate(VIS_cut = case_when(
    VIS_D0 > 15 ~ ">15",
    VIS_D0 <= 15 ~ "≤15"
  )) %>%
  mutate(VIS_cut = factor(VIS_cut, levels = c("≤15", ">15")))

# ADRB1 in lymphocytes
F_VIS_adrb1 <- df %>%
  select(ADR1_JS, VIS_cut) %>%
  filter(!is.na(ADR1_JS), !is.na(VIS_cut))

count_data_adrb1 <- F_VIS_adrb1 %>%
  group_by(VIS_cut) %>%
  summarise(n = n(), .groups = "drop")

p_value_adrb1 <- wilcox.test(log10(ADR1_JS) ~ VIS_cut, data = F_VIS_adrb1)$p.value
p_label_adrb1 <- paste("p =", format(p_value_adrb1, digits = 2, nsmall = 3))

plot_ADRB1 <- ggplot(F_VIS_adrb1, aes(x = VIS_cut, y = log10(ADR1_JS), fill = VIS_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_adrb1, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
            inherit.aes = FALSE, size = 4) +
  annotate("text", x = 1.5, y = max(log10(F_VIS_adrb1$ADR1_JS), na.rm = TRUE), label = p_label_adrb1, size = 5) +
  scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
  labs(title = "ADRB1 expression in lymphocytes by VIS score", x = "VIS score", y = "ADRB1 in lymphocytes at day 0 (log10 scale)") +
  theme_minimal() +
  theme(legend.position = "none")

# ADRB2 in lymphocytes
F_VIS_adrb2 <- df %>%
  select(ADR2_JS, VIS_cut) %>%
  filter(!is.na(ADR2_JS), !is.na(VIS_cut))

count_data_adrb2 <- F_VIS_adrb2 %>%
  group_by(VIS_cut) %>%
  summarise(n = n(), .groups = "drop")

p_value_adrb2 <- wilcox.test(log10(ADR2_JS) ~ VIS_cut, data = F_VIS_adrb2)$p.value
p_label_adrb2 <- paste("p =", format(p_value_adrb2, digits = 2, nsmall = 3))

plot_ADRB2 <- ggplot(F_VIS_adrb2, aes(x = VIS_cut, y = log10(ADR2_JS), fill = VIS_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_adrb2, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
            inherit.aes = FALSE, size = 4) +
  annotate("text", x = 1.5, y = max(log10(F_VIS_adrb2$ADR2_JS), na.rm = TRUE), label = p_label_adrb2, size = 5) +
  scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
  labs(title = "ADRB2 expression in lymphocytes by VIS score", x = "VIS score", y = "ADRB2 in lymphocytes at day 0 (log10 scale)") +
  theme_minimal() +
  theme(legend.position = "none")

# Combined plot per linfociti
plot_ADRB_lympho_combined <- plot_ADRB1 + plot_ADRB2 + 
  plot_layout(ncol = 2) + 
  plot_annotation(title = "ADRB1 and ADRB2 expression in lymphocytes stratified by VIS score")

print(plot_ADRB_lympho_combined)

ADRB ~ VIS Score at day 0 in Monocytes

# VIS_cut 
df <- df %>%
  mutate(VIS_cut = case_when(
    VIS_D0 > 15 ~ ">15",
    VIS_D0 <= 15 ~ "≤15"
  )) %>%
  mutate(VIS_cut = factor(VIS_cut, levels = c("≤15", ">15")))

### ADRB1 
F_VIS_adrb1_m <- df %>%
  select(ADR1m_JS, VIS_cut) %>%
  filter(!is.na(ADR1m_JS), !is.na(VIS_cut))

count_data_adr1m <- F_VIS_adrb1_m %>%
  group_by(VIS_cut) %>%
  summarise(n = n(), .groups = "drop")

p_value_adr1m <- wilcox.test(log10(ADR1m_JS) ~ VIS_cut, data = F_VIS_adrb1_m)$p.value
p_label_adr1m <- paste("p =", format(p_value_adr1m, digits = 2, nsmall = 3))

plot_ADRB1m <- ggplot(F_VIS_adrb1_m, aes(x = VIS_cut, y = log10(ADR1m_JS), fill = VIS_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_adr1m, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
            inherit.aes = FALSE, size = 4) +
  annotate("text", x = 1.5, y = max(log10(F_VIS_adrb1_m$ADR1m_JS), na.rm = TRUE), label = p_label_adr1m, size = 5) +
  scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
  labs(title = "ADRB1 expression in monocytes by VIS", x = "VIS score", y = "ADRB1 in monocytes at day 0 (log10 scale)") +
  theme_minimal() +
  theme(legend.position = "none")

### ADRB2 
F_VIS_adrb2_m <- df %>%
  select(ADR2m_JS, VIS_cut) %>%
  filter(!is.na(ADR2m_JS), !is.na(VIS_cut))

count_data_adr2m <- F_VIS_adrb2_m %>%
  group_by(VIS_cut) %>%
  summarise(n = n(), .groups = "drop")

p_value_adr2m <- wilcox.test(log10(ADR2m_JS) ~ VIS_cut, data = F_VIS_adrb2_m)$p.value
p_label_adr2m <- paste("p =", format(p_value_adr2m, digits = 2, nsmall = 3))

plot_ADRB2m <- ggplot(F_VIS_adrb2_m, aes(x = VIS_cut, y = log10(ADR2m_JS), fill = VIS_cut)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_adr2m, aes(x = VIS_cut, y = -2.5, label = paste0("n=", n)),
            inherit.aes = FALSE, size = 4) +
  annotate("text", x = 1.5, y = max(log10(F_VIS_adrb2_m$ADR2m_JS), na.rm = TRUE), label = p_label_adr2m, size = 5) +
  scale_fill_manual(values = c("≤15" = "#0072B2", ">15" = "#D55E00")) +
  labs(title = "ADRB2 expression in monocytes by VIS", x = "VIS score", y = "ADRB2 in monocytes at day 0 (log10 scale)") +
  theme_minimal() +
  theme(legend.position = "none")

# Combined
plot_ADRB_mono_combined <- plot_ADRB1m + plot_ADRB2m + 
  plot_layout(ncol = 2) + 
  plot_annotation(title = "ADRB1 and ADRB2 expression in monocytes stratified by VIS score")

# Combined plot per linfociti
plot_ADRB_lympho_combined <- plot_ADRB1 + plot_ADRB2 + 
  plot_layout(ncol = 2) + 
  plot_annotation(title = "ADRB1 and ADRB2 expression in lymphocytes stratified by VIS score")

print(plot_ADRB_mono_combined)

ADR ~ levosimendan

Figure 5: Association of Levosimendan and ADRB in Lymphocyte at day 3-5

#ADRB1

F5_lym <- df %>%
  dplyr::select(Levosimendan, ADR1_J3_J5) %>%
  filter(!is.na(Levosimendan), !is.na(ADR1_J3_J5)) %>%
  mutate(log_ADR1 = log10(ADR1_J3_J5))

# p-value
p_adr1 <- sprintf("%.4f", wilcox.test(log_ADR1 ~ Levosimendan, data = F5_lym, exact = FALSE)$p.value)

# count 
count_data_adr1 <- F5_lym %>%
  group_by(Levosimendan) %>%
  summarise(n = n(), .groups = "drop")

# Boxplot
Figure_5A <- ggplot(F5_lym, aes(x = Levosimendan, y = log_ADR1, fill = Levosimendan)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_adr1,
            aes(x = Levosimendan,
                y = min(F5_lym$log_ADR1, na.rm = TRUE) * 1.25,
                label = paste0("n=", n)),
            inherit.aes = FALSE,
            size = 4) +
  scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
  labs(title = "ADRB1 expression on T helper lymphocytes",
       x = "Levosimendan",
       y = "ADRB1 at D3–D5 of ECMO (log10 scale)") +
  annotate("text", x = 1.5, y = max(F5_lym$log_ADR1, na.rm = TRUE),
           label = paste("Mann–Whitney p =", p_adr1),
           size = 5, hjust = 0.5, color = "black") +
  theme_minimal() +
  theme(legend.position = "none")

#ADRB2

F5_lym_b <- df %>%
  dplyr::select(Levosimendan, ADR2_J3_J5) %>%
  filter(!is.na(Levosimendan), !is.na(ADR2_J3_J5)) %>%
  mutate(log_ADR2 = log10(ADR2_J3_J5))

# p-value
p_adr2 <- sprintf("%.4f", wilcox.test(log_ADR2 ~ Levosimendan, data = F5_lym_b, exact = FALSE)$p.value)

# count 
count_data_adr2 <- F5_lym_b %>%
  group_by(Levosimendan) %>%
  summarise(n = n(), .groups = "drop")

# Boxplot
Figure_5B <- ggplot(F5_lym_b, aes(x = Levosimendan, y = log_ADR2, fill = Levosimendan)) +
  geom_boxplot() +
  geom_point(position = position_jitter(width = 0.1), size = 2, alpha = 0.7) +
  geom_text(data = count_data_adr2,
            aes(x = Levosimendan,
                y = min(F5_lym_b$log_ADR2, na.rm = TRUE) * 1.25,
                label = paste0("n=", n)),
            inherit.aes = FALSE,
            size = 4) +
  scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
  labs(title = "ADRB2 expression on T helper lymphocytes",
       x = "Levosimendan",
       y = "ADRB2 at D3–D5 of ECMO (log10 scale)") +
  annotate("text", x = 1.5, y = max(F5_lym_b$log_ADR2, na.rm = TRUE),
           label = paste("Mann–Whitney p =", p_adr2),
           size = 5, hjust = 0.5, color = "black") +
  theme_minimal() +
  theme(legend.position = "none")

Figure_5A <- Figure_5A + theme(plot.title = element_blank())
Figure_5B <- Figure_5B + theme(plot.title = element_blank())

Figure5_tot_log <- Figure_5A + Figure_5B +
  plot_annotation(title = "Association between Levosimendan and ADRB expression (log scale) at D3–D5")

Figure5_tot_log

Figure 5m: Association of Levosimendan and ADRB in Monocyte at day of weaning

# ADRB1 
F5_mon_a <- df %>%
  dplyr::select(Levosimendan, ADR1m_J3_J5) %>%
  filter(!is.na(Levosimendan), !is.na(ADR1m_J3_J5))

# P-value
p_adr1m <- sprintf("%.4f", wilcox.test(ADR1m_J3_J5 ~ Levosimendan, data = F5_mon_a, exact = FALSE)$p.value)

# Count per group
count_data_adr1m <- F5_mon_a %>%
  group_by(Levosimendan) %>%
  summarise(n = n())

# Boxplot
Figure_5Am <- ggplot(F5_mon_a, aes(x = Levosimendan, y = log10(ADR1m_J3_J5), fill = Levosimendan)) +
  geom_boxplot() +
  geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
  geom_text(data = count_data_adr1m,
            aes(x = Levosimendan, y = min(log10(F5_mon_a$ADR1m_J3_J5), na.rm = TRUE) * 1.2,
                label = paste0("n=", n)),
            inherit.aes = FALSE,
            size = 4) +
  scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
  labs(title = "A",
       x = "Levosimendan",
       y = "Monocyte ADRB1 expression at D3–D5 (log10 scale)") +
  annotate("text", x = 1.5, y = max(log10(F5_mon_a$ADR1m_J3_J5), na.rm = TRUE),
           label = paste("Mann–Whitney p =", p_adr1m),
           size = 5, color = "black", hjust = 0.5) +
  theme_minimal() +
  theme(legend.position = "none")


# ADRB2 
F5_mon_b <- df %>%
  dplyr::select(Levosimendan, ADR2m_J3_J5) %>%
  filter(!is.na(Levosimendan), !is.na(ADR2m_J3_J5))

# P-value
p_adr2m <- sprintf("%.4f", wilcox.test(ADR2m_J3_J5 ~ Levosimendan, data = F5_mon_b, exact = FALSE)$p.value)

count_data_adr2m <- F5_mon_b %>%
  group_by(Levosimendan) %>%
  summarise(n = n())

# Boxplot
Figure_5Bm <- ggplot(F5_mon_b, aes(x = Levosimendan, y = log10(ADR2m_J3_J5), fill = Levosimendan)) +
  geom_boxplot() +
  geom_point(size = 2, alpha = 0.7, position = position_jitter(width = 0.1)) +
  geom_text(data = count_data_adr2m,
            aes(x = Levosimendan, y = min(log10(F5_mon_b$ADR2m_J3_J5), na.rm = TRUE) * 1.2,
                label = paste0("n=", n)),
            inherit.aes = FALSE,
            size = 4) +
  scale_fill_manual(values = c("NO" = "#56B4E9", "YES" = "#F0E442")) +
  labs(title = "B",
       x = "Levosimendan",
       y = "Monocyte ADRB2 expression at D3–D5 (log10 scale)") +
  annotate("text", x = 1.5, y = max(log10(F5_mon_b$ADR2m_J3_J5), na.rm = TRUE),
           label = paste("Mann–Whitney p =", p_adr2m),
           size = 5, color = "black", hjust = 0.5) +
  theme_minimal() +
  theme(legend.position = "none")

#Combined

Figure_5Am <- Figure_5Am + theme(plot.title = element_blank())
Figure_5Bm <- Figure_5Bm + theme(plot.title = element_blank())

Figure5_mon_tot <- Figure_5Am + Figure_5Bm +
  plot_annotation(title = "Association between Levosimendan and ADRB expression on monocytes at D3–D5 of ECMO",
                  tag_levels = "A")

Figure5_mon_tot

Results

ADRB2 expression at D3-D5 was significantly lower in patients who did not received Levosimendan compared to those who did. N difference where seen in monocytes.

Figure 6: Survival ~ ADR in Lymphocyte

# ADRB1
km_fit <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR1_D0_median, data = df)

Figure_6 <- ggsurvplot(km_fit, data = df, 
  pval = TRUE,               
  pval.coord = c(10, 0.1),   
  conf.int = FALSE,           
  risk.table = TRUE,         
  risk.table.height = 0.2,   
  risk.table.y.text = TRUE,  
  risk.table.title = "Number at Risk (number censored)",  
  xlim = c(0, 28),         
  break.time.by = 7,
  legend.title = " ",       
  palette = c("#E63946", "#457B9D"),
  ggtheme = theme_minimal(),
  legend.labs = c("High", "Low")
)

Figure_6 <- Figure_6 + ggtitle("28-day Survival Curves for ADRB1")

# ADRB2
km_fit <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR2_D0_median, data = df)

Figure_6.1 <- ggsurvplot(km_fit, data = df, 
  pval = TRUE,               
  pval.coord = c(10, 0.1),   
  conf.int = FALSE,           
  risk.table = TRUE,         
  risk.table.height = 0.2,   
  risk.table.y.text = TRUE,  
  risk.table.title = "Number at Risk (number censored)",  
  xlim = c(0, 28),         
  break.time.by = 7,
  legend.title = " ",       
  palette = c("#E63946", "#457B9D"),  
  ggtheme = theme_minimal(),
  legend.labs = c("High", "Low")
)

Figure_6.1 <- Figure_6.1 + ggtitle("28-day Survival Curves for ADRB2")

# Combine the plots (ADR1 and ADR2) side by side with the risk tables
Figure_6_combined <- plot_grid(Figure_6$plot, Figure_6$table, ncol = 1, align = "v", rel_heights = c(3, 1))
Figure_6.1_combined <- plot_grid(Figure_6.1$plot, Figure_6.1$table, ncol = 1, align = "v", rel_heights = c(3, 1))

# Remove titles from individual plots
Figure_6_combined <- Figure_6_combined + theme(plot.title = element_blank())
Figure_6.1_combined <- Figure_6.1_combined + theme(plot.title = element_blank())

# Combine both survival plots (ADR1 and ADR2) together with a title
Figure6_tot <- plot_grid(Figure_6_combined, Figure_6.1_combined, ncol = 2, align = "h")

# Add a title to the combined figure
Figure6_tot <- Figure6_tot + plot_annotation(title = "28-day Survival Curves for ADRB expressed in Lymphocytes")

Figure6_tot

Figure 6m: Survival ~ ADR in Monocyte

df_surv_m1 <- df %>%
  filter(!is.na(ADR1m_D0_median), !is.na(diff_days_J28), !is.na(Outcome_censored_28))

# ADRB1 on Monocytes
km_fit_m1 <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR1m_D0_median, data = df_surv_m1)


Fig_m1 <- ggsurvplot(km_fit_m1, data = df, 
  pval = TRUE,               
  pval.coord = c(10, 0.1),   
  conf.int = FALSE,           
  risk.table = TRUE,         
  risk.table.height = 0.2,   
  risk.table.y.text = TRUE,  
  risk.table.title = "Number at Risk (number censored)",  
  xlim = c(0, 28),         
  break.time.by = 7,
  legend.title = " ",       
  palette = c("#E63946", "#457B9D"),
  ggtheme = theme_minimal(),
  legend.labs = c("High", "Low")
)

Fig_m1 <- Fig_m1 + ggtitle("28-day Survival Curves for ADRB1")

# ADRB2 on Monocytes
km_fit_m2 <- survfit(Surv(diff_days_J28, Outcome_censored_28) ~ ADR2m_D0_median, data = df)

Fig_m2 <- ggsurvplot(km_fit_m2, data = df, 
  pval = TRUE,               
  pval.coord = c(10, 0.1),   
  conf.int = FALSE,           
  risk.table = TRUE,         
  risk.table.height = 0.2,   
  risk.table.y.text = TRUE,  
  risk.table.title = "Number at Risk (number censored)",  
  xlim = c(0, 28),         
  break.time.by = 7,
  legend.title = " ",       
  palette = c("#E63946", "#457B9D"),  
  ggtheme = theme_minimal(),
  legend.labs = c("High", "Low")
)

Fig_m2 <- Fig_m2 + ggtitle("28-day Survival Curves for ADRB2")

# Combine plots
Fig_m1_combined <- plot_grid(Fig_m1$plot, Fig_m1$table, ncol = 1, align = "v", rel_heights = c(3, 1))
Fig_m2_combined <- plot_grid(Fig_m2$plot, Fig_m2$table, ncol = 1, align = "v", rel_heights = c(3, 1))

Fig_m1_combined <- Fig_m1_combined + theme(plot.title = element_blank())
Fig_m2_combined <- Fig_m2_combined + theme(plot.title = element_blank())

Fig_monocyte_survival <- Fig_m1_combined + Fig_m2_combined +
  plot_annotation(title = "28-day Survival Curves for ADRB expressed in Monocytes")

Fig_monocyte_survival

Statistical analysis

To assess the prognostic significance of baseline β-adrenergic receptor expression, we stratified patients based on the median value of ADRB1 and ADRB2 measured at ECMO initiation (D0). Survival analysis was performed using Kaplan–Meier curves and the log-rank test to compare 28-day survival between patients with “High” vs “Low” expression. Patients were censored at the time of ECMO weaning, heart transplantation or LVAD implantation.

Results

Although the differences did not reach statistical significance, a trend was observed indicating that patients with higher ADRB expression in both lymphocytes and monocytes at ECMO initiation had lower 28-day survival probability. This pattern suggests a potential association between adrenergic receptor activation and worse short-term outcomes.

Citokines

plot_cytokine <- function(df, cytokine_name) {
  
  # Colomn name
  col_J0 <- paste0(cytokine_name, "_J0")
  col_J3_J5 <- paste0(cytokine_name, "_J3_J5")
  col_JS <- paste0(cytokine_name, "_JS")

  # 1. dataset
  data_plot <- df %>%
    select(ID, Outcome, all_of(c(col_J0, col_J3_J5, col_JS))) %>%
    mutate(
      Outcomes = factor(case_when(
        Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
        TRUE ~ as.character(Outcome)
      ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
    ) %>%
    pivot_longer(
      cols = c(all_of(c(col_J0, col_J3_J5, col_JS))),
      names_to = "cyto_time",
      values_to = "value"
    ) %>%
    mutate(
      time = case_when(
        cyto_time == col_J0 ~ 1,
        cyto_time == col_J3_J5 ~ 2,
        cyto_time == col_JS ~ 3,
        TRUE ~ NA_real_
      ),
      value_log = log10(value)
    )

  # 2. n
  count_data <- data_plot %>%
    filter(!is.na(value)) %>%
    group_by(cyto_time, Outcomes) %>%
    summarise(n = n(), .groups = "drop") %>%
    filter(!is.na(Outcomes))

  # 3. Models
  mod1 <- lmerTest::lmer(
    value_log ~ time * Outcomes + (1 + time | ID),
    data = data_plot
  )
  mod_no_interaction <- lmerTest::lmer(
    value_log ~ time + Outcomes + (1 + time | ID),
    data = data_plot
  )

  p_anova <- anova(mod_no_interaction, ddf = "Kenward-Roger")
  p_time <- format.pval(p_anova["time", "Pr(>F)"], digits = 3, eps = .001)
  p_outcome <- format.pval(p_anova["Outcomes", "Pr(>F)"], digits = 3, eps = .001)

  # 4. Plot
  p <- ggplot(data_plot, aes(x = cyto_time, y = value_log, fill = Outcomes)) +
    geom_boxplot(position = position_dodge(0.7), width = 0.6) +
    scale_y_log10() +
    labs(
      title = paste("Association between", cytokine_name, "and outcomes"),
      x = "Times of measurement",
      y = paste("Log10 of", cytokine_name, "(pg/mL)")
    ) +
    geom_text(data = count_data, aes(
      x = cyto_time,
      y = -min(data_plot$value_log, na.rm = TRUE) * 0.01,
      label = paste0("n=", n),
      group = Outcomes
    ), position = position_dodge(0.7), size = 4, vjust = 1) +
    scale_x_discrete(labels = setNames(
  c("implantation", "day 3 to 5", "explantation"),
  c(col_J0, col_J3_J5, col_JS)
)) +
    scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",
      "Bridge to transplant or LVAD" = "#377EB8",
      "Death" = "#F8766D"
    )) +
    annotate(
      "text",
      x = 1,
      y = max(data_plot$value_log, na.rm = TRUE) * 1.1,
      label = paste0("p[time] == ", p_time, "*';'*~p[outcomes] == ", p_outcome),
      parse = TRUE,
      size = 4,
      hjust = 0
    )

  print(p)
}

IL-6

plot_cytokine(df, "IL6")

IL-6 levels were lower in patients successfully weaned from ECMO, particularly at baseline (day 0) and at explantation. A significant effect of time was observed (p = 0.02), while the association with clinical outcomes showed a trend without reaching statistical significance (p = 0.055).

IL-33

plot_cytokine_linear <- function(df, cytokine_name = "IL-33") {
  col_J0 <- paste0(cytokine_name, "_J0")
  col_J3_J5 <- paste0(cytokine_name, "_J3_J5")
  col_JS <- paste0(cytokine_name, "_JS")

  data_plot <- df %>%
    select(ID, Outcome, all_of(c(col_J0, col_J3_J5, col_JS))) %>%
    mutate(
      Outcomes = factor(case_when(
        Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
        TRUE ~ as.character(Outcome)
      ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
    ) %>%
    pivot_longer(cols = all_of(c(col_J0, col_J3_J5, col_JS)),
                 names_to = "cyto_time", values_to = "value") %>%
    filter(!is.na(value)) %>%
    mutate(
      time = case_when(
        cyto_time == col_J0 ~ 1,
        cyto_time == col_J3_J5 ~ 2,
        cyto_time == col_JS ~ 3,
        TRUE ~ NA_real_
      )
    )

  count_data <- data_plot %>%
    group_by(cyto_time, Outcomes) %>%
    summarise(n = n(), .groups = "drop") %>%
    filter(!is.na(Outcomes))

  mod1 <- lmerTest::lmer(value ~ time * Outcomes + (1 + time | ID), data = data_plot)
  mod_no_interaction <- lmerTest::lmer(value ~ time + Outcomes + (1 + time | ID), data = data_plot)
  p_anova <- anova(mod_no_interaction, ddf = "Kenward-Roger")
  p_time <- format.pval(p_anova["time", "Pr(>F)"], digits = 3, eps = .001)
  p_outcome <- format.pval(p_anova["Outcomes", "Pr(>F)"], digits = 3, eps = .001)

  ggplot(data_plot, aes(x = cyto_time, y = value, fill = Outcomes)) +
    geom_boxplot(position = position_dodge(0.7), width = 0.6) +
    labs(
      title = "Association between IL-33 and outcomes",
      x = "Times of measurement",
      y = "IL-33 (pg/mL)"
    ) +
    geom_text(data = count_data, aes(
      x = cyto_time,
      y = min(data_plot$value, na.rm = TRUE) * 0.9,
      label = paste0("n=", n),
      group = Outcomes
    ), position = position_dodge(0.7), size = 4, vjust = 1) +
    scale_x_discrete(labels = setNames(
  c("implantation", "day 3 to 5", "explantation"),
  c(col_J0, col_J3_J5, col_JS)
)) +
    scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",
      "Bridge to transplant or LVAD" = "#377EB8",
      "Death" = "#F8766D"
    )) +
    annotate("text", x = 1, y = max(data_plot$value, na.rm = TRUE) * 1.1,
             label = paste0("p[time] == ", p_time, "*';'*~p[outcomes] == ", p_outcome),
             parse = TRUE, size = 4, hjust = 0)
}

plot_cytokine_linear(df, "IL-33")

Plasma concentrations of IL-33 were undetectable in the vast majority of patients, with only three samples showing quantifiable values above the detection threshold. Due to this distribution, IL-33 levels were visualized on a linear scale rather than log-transformed.

In addition to IL-33, we also measured IL-1β, IL-2 and IL-4; however, these cytokines were consistently undetectable across all samples, preventing any meaningful quantitative analysis. Similarly, plasma concentrations of IL-12p70, IFN-γ and Granzyme B were below the detection threshold in all patients. IL-17A was undetectable in the vast majority of samples, with only one patient showing a quantifiable value above the lower limit of detection. TNF-α was quantifiable in only two patients. Due to the pervasive non-detectability of these markers, no statistical comparison across timepoints or clinical outcomes was feasible.

IL-5

plot_cytokine(df, "IL5")

Among the cytokines measured, IL-5 concentrations appeared lower in patients successfully weaned from ECMO, although differences across timepoints and outcomes did not reach statistical significance.

IL-8

plot_cytokine(df, "IL8/CXCL8")

Also IL-8 concentrations appeared clear lower in patients successfully weaned from ECMO, although differences across timepoints and outcomes did not reach statistical significance.

IL-10

plot_cytokine <- function(df, cytokine_name) {
  
  # Column names for the three timepoints
  col_J0 <- paste0(cytokine_name, "_J0")
  col_J3_J5 <- paste0(cytokine_name, "_J3_J5")
  col_JS <- paste0(cytokine_name, "_JS")

  # 1. Prepare dataset
  data_plot <- df %>%
    select(ID, Outcome, all_of(c(col_J0, col_J3_J5, col_JS))) %>%
    mutate(
      Outcomes = factor(case_when(
        Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
        TRUE ~ as.character(Outcome)
      ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
    ) %>%
    pivot_longer(
      cols = c(all_of(c(col_J0, col_J3_J5, col_JS))),
      names_to = "cyto_time",
      values_to = "value"
    ) %>%
    mutate(
      time = case_when(
        cyto_time == col_J0 ~ 1,
        cyto_time == col_J3_J5 ~ 2,
        cyto_time == col_JS ~ 3,
        TRUE ~ NA_real_
      ),
      value_log = log10(value)
    )

  # 2. Count number per group
  count_data <- data_plot %>%
    filter(!is.na(value)) %>%
    group_by(cyto_time, Outcomes) %>%
    summarise(n = n(), .groups = "drop") %>%
    filter(!is.na(Outcomes))

  # 3. Fit models
  mod1 <- lmerTest::lmer(
    value_log ~ time * Outcomes + (1 + time | ID),
    data = data_plot
  )
  mod_no_interaction <- lmerTest::lmer(
    value_log ~ time + Outcomes + (1 + time | ID),
    data = data_plot
  )

  p_anova <- anova(mod_no_interaction, ddf = "Kenward-Roger")
  p_time_raw <- format.pval(p_anova["time", "Pr(>F)"], digits = 3, eps = .001)
  p_outcome_raw <- format.pval(p_anova["Outcomes", "Pr(>F)"], digits = 3, eps = .001)

  # Safe formatting for p-values with '<'
  p_time <- if (grepl("<", p_time_raw)) "\"<2e-16\"" else p_time_raw
  p_outcome <- if (grepl("<", p_outcome_raw)) "\"<2e-16\"" else p_outcome_raw

  # 4. Plot
  p <- ggplot(data_plot, aes(x = cyto_time, y = value_log, fill = Outcomes)) +
    geom_boxplot(position = position_dodge(0.7), width = 0.6) +
    scale_y_log10() +
    labs(
      title = paste("Association between", cytokine_name, "and outcomes"),
      x = "Times of measurement",
      y = paste("Log10 of", cytokine_name, "(pg/mL)")
    ) +
    geom_text(data = count_data, aes(
      x = cyto_time,
      y = -min(data_plot$value_log, na.rm = TRUE) * 0.01,
      label = paste0("n=", n),
      group = Outcomes
    ), position = position_dodge(0.7), size = 4, vjust = 1) +
    scale_x_discrete(labels = setNames(
      c("implantation", "day 3 to 5", "explantation"),
      c(col_J0, col_J3_J5, col_JS)
    )) +
    scale_fill_manual(values = c(
      "ECMO Weaning" = "#4DAF4A",
      "Bridge to transplant or LVAD" = "#377EB8",
      "Death" = "#F8766D"
    )) +
    annotate(
      "text",
      x = 1,
      y = max(data_plot$value_log, na.rm = TRUE) * 1.1,
      label = paste0("p[time] == ", p_time, "*';'*~p[outcomes] == ", p_outcome),
      parse = TRUE,
      size = 4,
      hjust = 0
    )

  print(p)
}

plot_cytokine(df, "IL10")

IL-10 levels showed a statistically significant change over time (p < 0.05), indicating dynamic modulation during the clinical course. However, no significant differences were observed across the defined outcome groups.

GDF-15

plot_cytokine(df, "GDF-15")

GDF-15 concentrations appeared lower in patients successfully weaned from ECMO, although differences across timepoints and outcomes did not reach statistical significance.

Results

Not all cytokines, whether pro-inflammatory or anti-inflammatory, were consistently measurable in our cohort. Based on the data, cytokine levels do not appear to be significantly associated with clinical outcomes. Interestingly, both pro-inflammatory and anti-inflammatory cytokines tend to decrease over time.

Citokynes Levels at Day 0 by Outcome Group

anti_infiamm <- c("IL10", "GDF-15")
pro_infiamm <- c("IL5", "IL6", "IL8/CXCL8")

plot_cytokines_group <- function(df, timepoint = "J3_J5") {

  cols_anti <- paste0(anti_infiamm, "_", timepoint)
  cols_pro <- paste0(pro_infiamm, "_", timepoint)
  
  df_sub <- df %>%
    select(ID, Outcome, all_of(c(cols_anti, cols_pro))) %>%
    mutate(
      Outcomes = factor(case_when(
        Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
        TRUE ~ Outcome
      ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
    )
  
  # Reshape long 
  df_long <- df_sub %>%
    pivot_longer(cols = -c(ID, Outcome, Outcomes),
                 names_to = "cytokine_time", values_to = "value") %>%
    mutate(
      cytokine = sub(paste0("_", timepoint), "", cytokine_time),
      # Negativ sign for anti-inflamm on left
      value_plot = ifelse(cytokine %in% anti_infiamm, -value, value),
      cytokine = factor(cytokine, levels = c(anti_infiamm, pro_infiamm))
    ) %>%
    filter(!is.na(value))
  
  # Colors
  colori <- c(
    "ECMO Weaning" = "#4DAF4A",
    "Bridge to transplant or LVAD" = "#377EB8",
    "Death" = "#F8766D"
  )
  
  # Plot
  ggplot(df_long, aes(x = value_plot, y = fct_rev(cytokine), fill = Outcomes)) +
    geom_bar(stat = "summary", fun = mean, position = "dodge", width = 0.7) +
    scale_fill_manual(values = colori) +
    scale_x_continuous(
      breaks = scales::pretty_breaks(n = 10),
      labels = abs,
      name = "Mean cytokine level (pg/mL)"
    ) +
    labs(title = paste("Cytokine levels at", timepoint, "by outcome group"),
         y = "Cytokine") +
    theme_minimal() +
    theme(
      legend.position = "right"
    )
}


plot_cytokines_group(df, "J0")

Citokynes Levels at Day 3-5 by Outcome Group

plot_cytokines_group(df, "J3_J5") 

Citokynes Levels at Day of Weaning by Outcome Group

plot_cytokines_group(df, "JS")

Citokynes Levels at Day 0 by Outcome Group, log scale

anti_infiamm <- c("IL10", "GDF-15")
pro_infiamm <- c("IL5", "IL6", "IL8/CXCL8")

plot_cytokines_group_log_vertical <- function(df, timepoint = "J3_J5") {
  
  cols_anti <- paste0(anti_infiamm, "_", timepoint)
  cols_pro <- paste0(pro_infiamm, "_", timepoint)
  
  df_sub <- df %>%
    select(ID, Outcome, all_of(c(cols_anti, cols_pro))) %>%
    mutate(
      Outcomes = factor(case_when(
        Outcome %in% c("Bridge to LVAD", "Bridge to Transplant") ~ "Bridge to transplant or LVAD",
        TRUE ~ Outcome
      ), levels = c("Death", "Bridge to transplant or LVAD", "ECMO Weaning"))
    )
  
  df_long <- df_sub %>%
    pivot_longer(cols = -c(ID, Outcome, Outcomes),
                 names_to = "cytokine_time", values_to = "value") %>%
    mutate(
      cytokine = sub(paste0("_", timepoint), "", cytokine_time),
      # Replace zeros or negatives before log10
      value_adj = ifelse(value <= 0, 1e-7, value),
      value_log = log10(value_adj),
      # Group label
      group = case_when(
        cytokine %in% anti_infiamm ~ "Anti-inflammatory",
        cytokine %in% pro_infiamm ~ "Pro-inflammatory",
        TRUE ~ "Other"
      ),
      # Custom sign flipping: IL10 positive, IL5 negative, other anti-inflam negative, pro positive
      value_plot = case_when(
        cytokine == "IL10" ~ value_log,
        cytokine == "IL5" ~ -value_log,
        group == "Anti-inflammatory" & cytokine != "IL10" ~ -value_log,
        TRUE ~ value_log
      ),
      cytokine = factor(cytokine, levels = c(anti_infiamm, pro_infiamm)),
      group = factor(group, levels = c("Anti-inflammatory", "Pro-inflammatory"))
    ) %>%
    filter(!is.na(value_log))
  
  colori <- c(
    "ECMO Weaning" = "#4DAF4A",
    "Bridge to transplant or LVAD" = "#377EB8",
    "Death" = "#F8766D"
  )
  
  # Plot base
  p <- ggplot(df_long, aes(x = value_plot, y = fct_rev(cytokine), fill = Outcomes)) +
    geom_bar(stat = "summary", fun = mean, position = position_dodge(width = 0.7), width = 0.6) +
    scale_fill_manual(values = colori) +
    scale_x_continuous(
      breaks = scales::pretty_breaks(n = 10),
      labels = function(x) round(abs(x), 2),
      name = "Mean log10 cytokine level (pg/mL)"
    ) +
    labs(title = paste("Log10 Cytokine levels at", timepoint, "by outcome group"),
         y = "Cytokine") +
    theme_minimal(base_size = 14) +
    theme(
      legend.position = "right",
      axis.text.y = element_text(size = 12),
      plot.margin = margin(t = 35, r = 20, b = 10, l = 20)
    )
  
  # Add top annotation for anti/pro groups with grid::annotation_custom
  library(grid)
  anti_label <- textGrob("Anti-inflammatory", x = unit(0.15, "npc"), y = unit(0.98, "npc"),
                        gp = gpar(fontsize = 14, fontface = "bold"))
  pro_label <- textGrob("Pro-inflammatory", x = unit(0.85, "npc"), y = unit(0.98, "npc"),
                       gp = gpar(fontsize = 14, fontface = "bold"))
  
  p + annotation_custom(anti_label) + annotation_custom(pro_label)
}

plot_cytokines_group_log_vertical(df, "J0")

At day 0, patients who later achieved ECMO weaning tended to have higher levels of anti-inflammatory cytokines such as IL-10. GDF-15 levels were comparable between patients who died and those who were weaned. In contrast, pro-inflammatory cytokines—particularly IL-6 and IL-8—were elevated in patients who died.

(For Antoine: For better graphical representation, I choose a logarithmic scale. Among all timepoints, I will focused primarily on day 0, aiming to identify potential predictors or markers that could help clinicians recognize patients with a more aggressive disease course and possibly guide management decisions, although currently no targeted therapy has been proven effective.)

# CK dosable
cyto_selected <- c("IL10", "GDF-15", "IL5", "IL6", "IL8/CXCL8")

anti_infiamm <- c("IL10", "GDF-15")

# At J0
cyto_cols <- paste0(cyto_selected, "_J0")

# Cutoff mediana + indosable
quantile_3groups_label_simple <- function(x, detection_limit = 0.001) {
  x_label <- ifelse(is.na(x) | x < detection_limit, "Indosabile", NA_character_)
  x_valid <- x[!(is.na(x) | x < detection_limit)]
  med_val <- median(x_valid, na.rm = TRUE)
  x_label[!(is.na(x) | x < detection_limit)] <- ifelse(x_valid <= med_val, "Low", "High")
  factor(x_label, levels = c("Indosabile", "Low", "High"))
}

# Groups High/Low ADRB1 e ADRB2 in lympho e mono (median)
df <- df %>%
  mutate(
    ADR1m_group = ifelse(ADR1m_J0 > median(ADR1m_J0, na.rm = TRUE), "High", "Low"),
    ADR2m_group = ifelse(ADR2m_J0 > median(ADR2m_J0, na.rm = TRUE), "High", "Low"),
    ADR1_group = ifelse(ADR1_J0 > median(ADR1_J0, na.rm = TRUE), "High", "Low"),
    ADR2_group = ifelse(ADR2_J0 > median(ADR2_J0, na.rm = TRUE), "High", "Low")
  )

df_quart <- df %>%
  select(ID, all_of(cyto_cols), ADR1m_group, ADR2m_group, ADR1_group, ADR2_group) %>%
  mutate(across(all_of(cyto_cols), ~ quantile_3groups_label_simple(.))) %>%
  filter(!is.na(ADR1m_group), !is.na(ADR2m_group), !is.na(ADR1_group), !is.na(ADR2_group))


make_long <- function(df, group_col, cell_type, adr_type) {
  df %>%
    select(ID, all_of(cyto_cols), all_of(group_col)) %>%
    pivot_longer(cols = all_of(cyto_cols), names_to = "cytokine", values_to = "level") %>%
    mutate(
      cytokine = sub("_J0$", "", cytokine),
      adr_group = .data[[group_col]],
      cell_type = cell_type,
      adr_type = adr_type
    ) %>%
    select(-all_of(group_col))
}

df_long_lymph_ADR1 <- make_long(df_quart, "ADR1_group", "Lymphocytes", "ADRB1")
df_long_lymph_ADR2 <- make_long(df_quart, "ADR2_group", "Lymphocytes", "ADRB2")
df_long_mono_ADR1 <- make_long(df_quart, "ADR1m_group", "Monocytes", "ADRB1")
df_long_mono_ADR2 <- make_long(df_quart, "ADR2m_group", "Monocytes", "ADRB2")

heatmap_data <- bind_rows(df_long_lymph_ADR1, df_long_lymph_ADR2, df_long_mono_ADR1, df_long_mono_ADR2) %>%
  mutate(
    cytokine = factor(cytokine, levels = cyto_selected),
    adr_group = factor(adr_group, levels = c("Low", "High")),
    cell_type = factor(cell_type, levels = c("Lymphocytes", "Monocytes")),
    adr_type = factor(adr_type, levels = c("ADRB1", "ADRB2"))
  )

palette_colors <- c(
  "Indosabile" = "white",
  "Low" = "forestgreen",
  "High" = "firebrick"
)


# Plot 
ggplot(heatmap_data, aes(x = adr_group, y = cytokine, fill = quartile)) +
  # Heatmap dati
  geom_tile(data = heatmap_data, aes(x = adr_group, y = cytokine, fill = level), color = "grey80", width = 0.8) +
  facet_grid(cell_type ~ adr_type) +
  scale_fill_manual(values = palette_colors, na.value = "grey90") +
  scale_x_discrete(expand = expansion(add = c(1, 0.1))) +
  coord_cartesian(clip = "off") +
  labs(
    title = "Heatmap of Selected Cytokine Levels at Day 0 by ADRB Expression Groups",
    x = "ADRB Expression Group",
    y = "Cytokine",
    fill = "Level"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(size = 12),
    axis.text.y = element_text(color = "black", angle = 45, hjust = 1),
    panel.grid = element_blank(),
    strip.background = element_blank(),
    strip.text.x = element_text(hjust = 0.5),
    strip.text.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)
  )

# I don't understand why IL5 and IL8 are white